The Interstellar Medium in High-redshift Submillimeter Galaxies as Probed by Infrared Spectroscopy∗

Julie L. Wardlow; Asantha Cooray; Willow Osage; Nathan Bourne; David Clements; Helmut Dannerbauer; Loretta Dunne; Simon Dye; Steve Eales; Duncan Farrah; Cristina Furlanetto; Edo Ibar; Rob Ivison; Steve Maddox; Michał M. Michałowski; Dominik Riechers; Dimitra Rigopoulou; Douglas Scott; Matthew W. L. Smith; Lingyu Wang; Paul van der Werf; Elisabetta Valiante; Ivan Valtchanov; Aprajita Verma

doi:10.3847/1538-4357/837/1/12

1. Introduction

Submillimeter galaxies (SMGs), selected from their high flux densities at submillimeter wavelengths, are the highest luminosity dusty star-forming galaxies and have redshift distributions peaking at $z\simeq 2$ with a tail out to $z\simeq 6$ (e.g., Chapman et al. 2005; Wardlow et al. 2011; Riechers et al. 2013; Dowell et al. 2014; Simpson et al. 2014; Asboth et al. 2016). They have intrinsic far-infrared (IR) luminosities ≳10¹² ${L}_{\odot }$ , equivalent to local ultraluminous infrared galaxies (ULIRGs), and are typically dominated by star formation rather than AGN emission (e.g., Alexander et al. 2005; Valiante et al. 2007; Pope et al. 2008; Menéndez-Delmestre et al. 2009; Laird et al. 2010; Wang et al. 2013). SMGs with fluxes down to ∼1 mJy at 850 μm contribute up to 20% of the cosmic star formation rate density at z = 2 (e.g., Wardlow et al. 2011; Swinbank et al. 2014). See Blain et al. (2002) and Casey et al. (2014) for reviews.

The extreme star formation rates of SMGs (up to ∼1000 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ ) and their gas-depletion times suggest that their star formation is episodic and that they are observed in a short-lived (timescales ∼100 Myr) burst phase (e.g., Bothwell et al. 2013). Both mergers and secular processes have been invoked as the triggers of these starbursts (e.g., Elbaz et al. 2011; Alaghband-Zadeh et al. 2012; Hayward et al. 2013; Menéndez-Delmestre et al. 2013; Cowley et al. 2015; Narayanan et al. 2015). Given the limited data, discourse is ongoing. A related issue is whether the star formation in SMGs proceeds like that in local ULIRGs (e.g., Daddi et al. 2010; Genzel et al. 2010), or whether the so-called "mode" of star formation proceeds more similarly to local sub-LIRGs or quiescently star-forming galaxies (e.g., Farrah et al. 2008; Pope et al. 2008; Elbaz et al. 2011; Krumholz et al. 2012), where it is typically extended over larger regions. The majority of local ULIRGs occur in interacting or merging systems (e.g., Sanders & Mirabel 1996; Farrah et al. 2001; Veilleux et al. 2002) but hints are beginning to emerge that SMGs may have a lower merger fraction (e.g., Tacconi et al. 2008, 2010; Rodighiero et al. 2011). There is also some evidence that the star-forming regions in SMGs may be more spatially extended than in local ULIRGs, suggestive of star formation proceeding in a sub-LIRG mode (e.g., Tacconi et al. 2006; Younger et al. 2008; Swinbank et al. 2010; Ivison et al. 2011; Riechers et al. 2011b; Ikarashi et al. 2015; Simpson et al. 2015), although recent lensing studies tend to measure smaller sizes than unlensed results. (e.g., Bussmann et al. 2013; Calanog et al. 2014).

Different star formation triggers, modes, and AGN contributions impact the ISM of galaxies and consequently manifest in the relative strengths of fine structure emission lines. Thus, observations of fine structure lines are crucial to investigate these aspects of SMGs. However, the dust that drives their extreme far-IR luminosities also makes observations at optical and near-IR wavelengths challenging, and renders standard excitation tracers inaccessible. Indeed, mid- and far-IR spectroscopy is the only way to probe the ISM in the inner, most highly extincted regions ( ${A}_{V}\gtrsim 6$ –10 mag). The limited wavelength coverage and sensitivity of previous mid-IR spectrographs (e.g., Spitzer/IRS, ISO/SWS, ISO/LWS) precluded observations of mid-IR fine structure lines for high-redshift galaxies prior to Herschel. Even with the enhanced sensitivity of Herschel, observations are limited to the brightest galaxies—primarily gravitationally lensed SMGs. Indeed, to date, only a handful of observations of the [O iv]26 μm, [S iii]33 μm, [Si ii]34 μm, [O iii]52 μm, [N iii]57 μm, or [O i]63 μm IR fine-structure lines have been observed in high-redshift galaxies, the majority taken with Herschel (Ivison et al. 2010a; Sturm et al. 2010; Valtchanov et al. 2011; Coppin et al. 2012; Bothwell et al. 2013; Brisbin et al. 2015, see also Carilli & Walter 2013 for a review of gas tracers in high-redshift galaxies).

In this paper, we present Herschel/PACS (Pilbratt et al. 2010; Poglitsch et al. 2010) observations of the [O iv]26 μm, [S iii]33 μm, [Si ii]34 μm, [O iii]52 μm, [N iii]57 μm, and [O i]63 μm fine structure transitions, and the molecular hydrogen rotational lines ${{\rm{H}}}_{2}$ S(0) (28 μm) and ${{\rm{H}}}_{2}$ S(1) (17 μm), in 13 strongly gravitationally lensed SMGs at redshifts 1.03–3.27 targeted by our Herschel Open Time program. These emission lines were selected to probe a range of ISM conditions, in terms of ionization potential and critical density. They correspond to different excitation mechanisms in photo-dominated regions (PDRs), H ii regions, shocks, and X-ray dominated regions (XDRs). We supplement these data with archival observations of the same IR emission lines from a further 32 SMGs (lensed and unlensed) at z = 1.1–4.2 from eight additional PACS observing programs. To complement the spectroscopy, we also obtained Herschel-PACS 70 and 160 μm photometry of the 13 original targets, which supplements the existing far-IR photometry of these lensed SMGs and is used to improve the SED fits. After the warm-up of Herschel, such spectroscopy will not again be attainable at high redshifts until the launch of facilities such as SPICA, FIRSPEX, or the Far-Infrared Surveyor. Thus, this paper represents one of the few studies of rest-frame mid-IR spectroscopy at high-redshifts in the present era, and provides important data for the planning of the observing strategies for these future missions.

In Section 2, we describe the observations and data reduction. Section 3, contains the analysis and discussion, including SED fits, emission line measurements, and ISM modeling. Finally, our conclusions are presented in Section 4. Throughout this paper, we use ΛCDM cosmology with ${{\rm{\Omega }}}_{{\rm{M}}}=0.27$ , ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$ and ${H}_{0}=71\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ .

2. Observations and Data Reduction

In this paper, we first analyze PACS observations of sources targeted by our Herschel program, OT2_jwardlow_1, as described in Section 2.1. We later combine these with archival spectroscopy for additional SMGs, as described in Section 2.4.

2.1. Targeted Sample Selection

The parent sample of the 13 galaxies targeted by OT2_jwardlow_1 for PACS photometry and spectroscopy are candidate strongly gravitationally lensed galaxies identified in the Herschel H-ATLAS (Eales et al. 2010) and HerMES (Oliver et al. 2012) surveys due to their brightness at 500 μm ( ${S}_{500}\geqslant 100$ mJy; e.g., Negrello et al. 2010; Wardlow et al. 2013; Nayyeri et al. 2016; Negrello et al. 2017). Extensive follow-up programs, including CO spectroscopy (e.g., Frayer et al. 2011; Harris et al. 2012; D. Riechers et al. 2017, in preparation), high-resolution (sub)millimeter and radio interferometry (e.g., Bussmann et al. 2013), high-resolution near-IR imaging (e.g., Wardlow et al. 2013; Calanog et al. 2014; Negrello et al. 2014), deep optical, near- and mid-IR photometry (e.g., Fu et al. 2013), and spectroscopy (e.g., Wardlow et al. 2013), are supplementing the ancillary data coverage of many of these systems.

The subset of gravitationally lensed Herschel-selected galaxies that are targeted here are presented in Table 1. The targeted galaxies were selected to have confirmed (multiple-line) CO spectroscopic redshifts, as well as ${S}_{250}\geqslant 100$ mJy and 70 μm fluxes predicted to be $\geqslant 5$ mJy, based on fitting Arp 220 and M 82 SEDs (Silva et al. 1998) to the available long wavelength data. The latter two requirements were motivated by the sensitivity of PACS, and the former is necessary to tune the spectroscopic observations (although many of the redshifts are from broadband instruments used for line searches, which can have up to $\sim 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ spectral resolution). PACS spectroscopy of six additional Herschel H-ATLAS and HerMES gravitationally lensed galaxies, as well as other high-redshift galaxies, were observed in a separate program and will be presented in A. Verma et al. (2017, in preparation). They are also included here in our stacking analyses (see Section 2.4).

Table 1. Positions, Redshifts, and the Lensing Amplifications of the Target Galaxies

Target	Short Names	${z}_{\mathrm{source}}$	${z}_{\mathrm{lens}}$	Magnification^a	References^b	OBSIDs^c
H-ATLAS J142935.3−002836	G15v2.19,	1.027	0.218	9.7 ± 0.7^a	C14, M14, N16	134225916[2, 3], 134226146[8, 9],
	G15.DR1.14					1342248369
H-ATLAS J085358.9+015537	G09v1.40,	2.089	...	15.3 ± 3.5	B13, C14, S16, Y16,	134225565[2, 3], 134225495[3–6],
	G09.DR1.35				N16	1342254283
H-ATLAS J115820.2−013753	G12v2.257,	2.191	...	13.0 ± 7.0	H12, N16	13422580[78–81], 134225725[1, 2],
	G12.DR1.379					1342257277
H-ATLAS J133649.9+291801	NGP.NA.144	2.202	...	4.4 ± 0.8	O13, H12, B13, N16	134225932[4, 5], 134225728[3–8]
H-ATLAS J134429.4+303036	NGP.NA.56	2.302	0.672	11.7 ± 0.9	H12, B13, Y16, N16	134225932[8, 9], 134225961[2–5],
						134225779[7, 8], 1342257289
1HerMES S250 J022016.5−060143	HXMM01	2.307	0.654	1.5 ± 0.3	B13, F13, W13,	134226195[7, 8], 1342262548,
					C14, B15	13422626[59, 60], 1342262769,
						1342263495
H-ATLAS J084933.4+021443	G09v1.124,	2.410	0.348	2.8 ± 0.2	H12, B13, C14, I14,	134225473[5–8], 13422549[57–60]
	G09.DR1.131				Y16, N16	1342254283
H-ATLAS J141351.9−000026	G15v2.235,	2.479	0.547	1.8 ± 0.3	H12, B14, C14, N16	13422591[58–61], 134226147[1, 2],
	G15.DR1.265					1342262532, 1342262041
H-ATLAS J091840.8+023047	G09v1.326,	2.581	...	1	H12, B13, C14, N16	134225564[6–9], 1342254933,
	G09.DR1.437					1342255740
H-ATLAS J133008.4+245900	NGP.NB.78	3.111	0.428	13.0 ± 1.5	O13, B13, C14, Y16,	134225932[0–3], 134225728[0–2]
					N16, Rp
H-ATLAS J113526.3−014605	G12v2.43,	3.128	...	2.8 ± 0.4	GY05, B13, C14,	13422571[09–12], 134225724[5–7],
	G12.DR1.80				Y16, N16	1342256482
H-ATLAS J114637.9−001132	G12v2.30,	3.259	1.225	9.5 ± 0.6	O13, F12, H12,	134225710[1–4], 13422572[48–50],
	G12.DR1.33				B13, C14, N16	1342256949, 1342257276
1HerMES S250 J143330.8+345439	HBoötes01	3.274	0.590	4.5 ± 0.4	B13, W13, C14, Rp	134225952[1–4], 13422620[39, 40],
						1342257689

Notes.

^aThe magnifications used are, with the exception of G15v2.19, for the far-IR continuum and measured from high-resolution submillimeter data (mostly observed-frame 850 μm with the SMA or ALMA). Section 3.7.1 includes further discussion of the effects of differential magnification. For G15v2.19, we use the magnification of the CO(4-3) line, which is the data closest in wavelength to our observations with lens modeling. ^bB13: Bussmann et al. (2013), B15: Bussmann et al. (2015), C14: Calanog et al. (2014), F12: Fu et al. (2012), GY05: Gladders & Yee (2005), H12: Harris et al. (2012), I13: Ivison et al. (2013), O13: Omont et al. (2013), M14: Messias et al. (2014), N16: Negrello et al. (2017), Rp: D. Riechers et al. (2017, in preparation), S16: Serjeant (2016), W13: Wardlow et al. (2013), Y16: Yang et al. (2016). ^cOBSID are the Herschel observation identification number(s) for the program OT2_jwardlow_1, used to identify the photometric and spectroscopic observation of each target in the Herschel archive.

Download table as: ASCII Typeset image

Table 2. PACS 70 and 160 μm Photometry and Derived Far-IR Fluxes and Luminosities

Name	S₇₀^a	S₁₆₀	${L}_{\mathrm{FIR}}$ ^b	${FIR}$ ^b
	(mJy)	(mJy)	( ${10}^{13}\,{L}_{\odot }$ )	( ${10}^{-15}\,{\mathrm{Wm}}^{-2}$ )
G15v2.19	316 ± 16	1077 ± 55	3.2 ± 0.3	16.6 ± 0.7
G09v1.40	$\lt 9$	279 ± 16	4.6 ± 0.3	4.0 ± 0.1
G12v2.257	15 ± 4	147 ± 11	1.7 ± 0.3	1.3 ± 0.1
GP.NA.144	11 ± 3	177 ± 12	${4.0}_{-0.3}^{+0.4}$	3.1 ± 0.1
NGP.NA.56	14 ± 3	303 ± 18	${7.2}_{-0.3}^{+0.4}$	5.0 ± 0.1
HXMM01	10 ± 3	123 ± 12	2.9 ± 0.3	2.0 ± 0.1
G09v1.124	16 ± 4	169 ± 11	4.4 ± 0.3	2.7 ± 0.1
G15v2.235	$\lt 11$	115 ± 9	3.6 ± 0.3	2.1 ± 0.1
G09v1.326	$\lt 8$	126 ± 9	3.1 ± 0.4	1.6 ± 0.1
NGP.NB.78	40 ± 3	210 ± 12	8.2 ± 0.7	2.7 ± 0.1
G12v2.43	16 ± 3	219 ± 12	${8.9}_{-0.6}^{+0.7}$	2.8 ± 0.1
G12v2.30	30 ± 4	228 ± 13	${11.2}_{-0.6}^{+0.7}$	3.3 ± 0.1
HBoötes01	$\lt 4$	81 ± 6	5.6 ± 0.5	1.7 ± 0.1

Notes. All measurements are apparent values (i.e., no corrections have been made for the lensing amplification).

^a $3\sigma$ $3\sigma$ upper limits are presented for undetected sources. ^b ${L}_{\mathrm{FIR}}$ ${L}_{\mathrm{FIR}}$ and ${FIR}$

${FIR}$ , measured from the modified blackbody fits in Figure 1 are the far-IR luminosity (40–500 μm) and continuum flux (42.5–122.5 μm), respectively (Section 3.1).

Download table as: ASCII Typeset image

2.2. Herschel-PACS Spectroscopy

The emission lines that were targeted vary from galaxy to galaxy, due to the redshift range of the sources and the PACS spectral coverage and sensitivity. In this section, we discuss the observations of the targeted sample of Herschel lensed SMGs (the data processing is the same for the archival data; Section 2.4). All of the targeted galaxies (Section 2.1) had between three and eight lines observed, with a median of five, from the [O iv]26, [S iii]33, [Si ii]34, [O iii]52, [N iii]57, and [O i]63 fine-structure transitions, and the molecular hydrogen rotational lines ${{\rm{H}}}_{2}$ S(0) and ${{\rm{H}}}_{2}$ S(1). The [Fe ii]26 transition is serendipitously included in the wavelength coverage of the [O iv]26 observations. The breakdown of the lines that were observed for each galaxy is shown in Table 3.

Table 3. Spectral Measurements

	[O i]63 μm			[S iii]33 μm			[Si ii]34 μm
Name	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b
	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)
G15v2.19	$\lt 35.3$	128.07	1044 ± 157	⋯	⋯	⋯	⋯	⋯	⋯
G09v1.40	$\lt 81.0$	195.29	$\lt 549$	⋯	⋯	⋯	$\lt 21.2$	107.61	$\lt 79$
G12v2.257	⋯	⋯	⋯	$\lt 10.5$	106.84	$\lt 39$	$\lt 10.9$	111.09	$\lt 42$
NGP.NA.144	$\lt 44.4$	202.30	$\lt 312$	$\lt 20.6$	107.21	$\lt 76$	$\lt 19.8$	111.47	$\lt 76$
NGP.NA.56	$\lt 34.0$	208.62	$\lt 245$	$\lt 24.5$	110.55	$\lt 94$	$\lt 27.0$	114.96	$\lt 107$
HXMM01	$\lt 12.7$	208.94	$\lt 92$	$\lt 13.7$	110.72	$\lt 52$	$\lt 15.6$	115.13	$\lt 62$
G09v1.124	⋯	⋯	⋯	$\lt 15.8$	114.17	$\lt 62$	$\lt 17.9$	118.72	$\lt 73$
G15v2.235	⋯	⋯	⋯	$\lt 12.8$	116.48	$\lt 51$	$\lt 11.4$	121.12	$\lt 47$
G09v1.326	⋯	⋯	⋯	$\lt 13.4$	119.96	$\lt 55$	⋯	⋯	⋯
NGP.NB.78	⋯	⋯	⋯	$\lt 19.1$	137.64	128 ± 91	$\lt 17.1$	143.12	151 ± 84
G12v2.43	⋯	⋯	⋯	$\lt 17.3$	138.21	103 ± 82	$\lt 15.3$	143.71	119 ± 76
G12v2.30	⋯	⋯	⋯	$\lt 32.4$	142.60	$\lt 160$	$\lt 29.7$	148.27	$\lt 152$
HBootes01	⋯	⋯	⋯	$\lt 11.2$	143.10	$\lt 55$	$\lt 12.3$	148.80	$\lt 63$
Mean Stack^c	1.0 ± 0.3	63.16	...	<0.9	⋯	...	2.8 ± 0.4	34.83	...

	[O iii]52 μm			[N iii]57 μm

Name	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b
	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)

G15v2.19	⋯	⋯	⋯	$\lt 32.9$	116.19	913 ± 132
G09v1.40	$\lt 45.4$	160.14	$\lt 252$	$\lt 18.6$	177.18	220 ± 114
G12v2.257	$\lt 14.3$	165.33	$\lt 82$	⋯	⋯	⋯
NGP.NA.144	10.6 ± 3.2	165.86	$\lt 134$	$\lt 22.8$	183.54	$\lt 145$
NGP.NA.56	$\lt 23.2$	171.08	241 ± 137	$\lt 28.5$	189.27	342 ± 187
HXMM01	$\lt 7.2$	171.34	129 ± 42	⋯	⋯	⋯
G09v1.124	$\lt 27.3$	176.67	$\lt 167$	⋯	⋯	⋯
G15v2.235	$\lt 24.5$	180.25	$\lt 153$	⋯	⋯	⋯
G09v1.326	$\lt 23.0$	185.64	$\lt 148$	⋯	⋯	⋯
NGP.NB.78	⋯	⋯	⋯	⋯	⋯	⋯
G12v2.43	⋯	⋯	⋯	⋯	⋯	⋯
G12v2.30	⋯	⋯	⋯	⋯	⋯	⋯
HBootes01	⋯	⋯	⋯	⋯	⋯	⋯
Mean Stack^c	<0.9	⋯	...	1.9 ± 0.6	57.38	...

	[O iv]26 μm^d and [Fe ii]26 μm^d

Name	[O iv] flux	[O iv] ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b	[Fe ii] flux	[Fe ii] ${\lambda }_{\mathrm{obs}}$ ^a
	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)	(10⁻¹⁸ Wm⁻²)	( μm)

G15v2.19	⋯	⋯	⋯	⋯	⋯
G09v1.40	<90.8	80.03	$\lt 252$	<91.2	80.33
G12v2.257	<51.4	82.61	$\lt 147$	<51.6	82.93
NGP.NA.144	<66.2	82.90	226 ± 190	<66.5	83.21
NGP.NA.56	<138.8	85.49	$\lt 411$	<139.3	85.81
HXMM01	7.7 ± 2.5	85.70	$\lt 66$	<22.6	85.94
G09v1.124	<95.9	88.28	$\lt 293$	<96.2	88.62
G15v2.235	<56.0	90.07	$\lt 175$	<56.2	90.41
G09v1.326	<62.6	92.76	$\lt 201$	<62.8	93.12
NGP.NB.78	⋯	⋯	⋯	⋯	⋯
G12v2.43	<7.0	106.87	53 ± 26	<7.0	107.28
G12v2.30	<10.7	110.27	56 ± 40	<10.7	110.68
HBootes01	⋯	⋯	⋯	⋯	⋯
Mean Stack^c	<4.0	⋯	⋯	<4.0	⋯


	${{\rm{H}}}_{2}$ S(0)			${{\rm{H}}}_{2}$ S(1)

Name	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b	Line flux	Line ${\lambda }_{\mathrm{obs}}$ ^a	Continuum^b
	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)	(10⁻¹⁸ Wm⁻²)	( μm)	(mJy)

G15v2.19	$\lt 94.0$	57.20	$\lt 186$	⋯	⋯	⋯
G09v1.40	<31.7	87.22	$\lt 96$	⋯	⋯	⋯
G12v2.257	⋯	⋯	⋯	⋯	⋯	⋯
NGP.NA.144	⋯	⋯	⋯	$\lt 69.9$	54.55	<132
NGP.NA.56	<41.5	93.18	$\lt 134$	$\lt 83.2$	56.25	<162
HXMM01	⋯	⋯	⋯	$\lt 44.4$	56.33	<86
G09v1.124	⋯	⋯	⋯	$\lt 54.2$	58.09	<109
G15v2.235	⋯	⋯	⋯	$\lt 31.1$	59.26	<64
G09v1.326	⋯	⋯	⋯	⋯	⋯	⋯
NGP.NB.78	<10.8	116.01	107 ± 43	⋯	⋯	⋯
G12v2.43	<9.9	116.49	86 ± 40	⋯	⋯	⋯
G12v2.30	<16.0	120.18	$\lt 66$	$\lt 34.2$	72.55	<86
HBootes01	<6.2	120.61	$\lt 26$	⋯	⋯	⋯
Mean Stack^c	<3.8	⋯	⋯	<4.2	⋯	⋯

Notes. We give $3\sigma$ upper limits for lines that are not detected above the $3\sigma$ significance level. Parameters for lines without observations are left blank.

^aFor detected lines, the wavelength corresponds to the measured (observed frame) position of the line. Otherwise, the expected (observed frame) wavelength is given, based on the nominal redshifts in Table 1. ^bThe continuum flux measured adjacent to the emission line. ^cThe bolded values are the stack values and are discussed in Section 3.3. ^dThe [O iv]26 and [Fe ii]26 lines occur close together in a single spectrum and are therefore fit simultaneously.

Download table as: ASCIITypeset images: 1 2

The data were taken in "range scan" mode with small chop/nod throws for background subtraction. With the exceptions of G12v2.30 and G12v2.43, the [O iv]26 lines were observed in the second order of the [O iii]52 observations. For G12v2.30 and G12v2.43, the [O iii]52 line is redshifted beyond the PACS wavelength range. In those cases, [O iv]26 was observed separately.

The data were reduced using the Herschel Interactive Processing Environment (Ott et al. 2010; hipe) v12.1.0 with version 65.0 of the PACS calibration tree.²¹ Data processing is based on the hipe v12.1.0 ipipe Background Normalization data reduction script for "chop/nod range scan" data. This procedure is optimized for faint sources and uses the off-source positions to perform the background subtraction and calibrate the detector response. During flat fielding, we set the "upsample factor" to 1 (and use the default "oversample" of 2) to avoid introducing correlated noise, and mask the wavelength regions where spectral lines are expected. The final spectra are binned to be Nyquist sampled at the native PACS resolution, and are shown in Appendix B. For the targets that are marginally resolved in the PACS photometry²² (Section 2.3), we applied the hipe extended source correction (assuming sizes measured at 70 μm); otherwise, we applied the standard point source correction during the extraction of the 1D spectra.

PACS always takes second-order spectroscopy, which, with the exception of the [O iv]26 and [O iii]52 observations described above, are not expected to include any additional transition lines. This is because no bright transitions of the background SMGs lie in the second-order wavelength ranges, and the foreground lensing galaxies are IR faint. Nevertheless the second-order data were reduced and extracted following the same procedure. As anticipated, no additional transitions were found. The continuum measurements (or limits) from these spectra are not deep enough to provide additional robust constraints on the SEDs. Therefore, the second-order data (with the exception of the paired [O iii]52 and [O iv]26 observations) are excluded from further examination.

2.3. Herschel-PACS Photometry

To supplement the spectroscopy, we also obtained simultaneous 70 and 160 μm mini-scan maps of each of the target lensed SMGs. Observations were taken at the nominal scan speed of 20''/s, with 3' scan legs, separated by 4'' cross-scan steps. For photometric fidelity, at least two orthogonal scans of each source were made. For the fainter targets, additional scan pairs were obtained to increase the observation depths.

The data were processed from level 0 using hipe v12.1.0 with version 65.0 of the PACS calibration tree. We employed standard Herschel data-reduction procedures, utilizing the standard ipipe script for scan maps containing point or marginally extended sources. The cross-scans were combined during reduction, and we iteratively filtered using a signal-to-noise (S/N) threshold to mask the sources during filtering. The final maps are each $\sim 3\buildrel{\,\prime}\over{.} 5\times 7\buildrel{\,\prime}\over{.} 5$ in size with coverage $\geqslant 90 \%$ of the maximum in the central $\sim 0\buildrel{\,\prime}\over{.} 5\times 1^{\prime}$ area.

PACS photometry is measured in 18 apertures, with radii from 2 to 50'', using the "annularSkyAperturePhotometry" task within hipe. Each measurement is corrected for the encircled energy fractions using PACS responsivity version 7. The uncertainties in the flux density measurements are determined from the dispersion in 1000 samples of the total flux in the same number of randomly selected pixels as included in each aperture, with pixels containing sources or those with $\lt 50 \%$ of maximum coverage excluded from selection. We then determine the "total" flux density and uncertainty for each target by fitting a curve of growth to the aperture fluxes and adding 5% calibration uncertainty.²³ These total flux densities are presented in Table 2, where we also include PACS 100 μm data from R. George et al. (2017, in preparation; Herschel program OT1_rivison_1; see also George 2015) and the publications presented in Table 1. George (2015) also includes 160 μm data from OT1_rivison_1, although their flux measurements can be 1–2σ lower than those presented here, because point sources are assumed. For HXMM01, the PACS 70 and 160 μm photometry was independently reduced and measured in Fu et al. (2013). Our measurements are consistent with those results, and we include the Fu et al. (2013) 100 μm photometry in the SED fits (Section 3.1) and Table 2. HBoötes01 has PACS 100 μm data from HerMES GTO time, which are also included here.

2.4. Archival Sample and Data

To identify additional SMGs with IR spectroscopy, we searched the successful Herschel proposals²⁴ for those targeting high-redshift star-forming galaxies (i.e., excluding AGN and QSOs) for PACS spectroscopy. Having identified likely programs, we next searched the Herschel Science Archive for those with SMGs as targets and retained the observations of IR emission lines that overlap with those studied by our own program (Section 2.2). This search resulted in spectroscopy for an additional 32 SMGs at CO or optical spectroscopic redshifts of 1.1 to 4.2. Most of these additional SMGs are gravitationally lensed because the PACS sensitivity means that only the apparently brightest sources can be observed. These archival observations covered between one and seven emission lines per galaxy. The full list of archival targets and data included in our analyses are presented in Table 5. The archival sources are broadly consistent with the main SMG population and the individually targeted galaxies, in terms of the IR-luminosity and redshift distributions, with IR emission being dominated by star formation. This archival sample includes LESS SMGs (Coppin et al. 2012), lensed HerMES and H-ATLAS sources from a similar followup program to this (OT1_averma_1; A. Verma et al. 2017, in preparation), lensed SPT sources (Vieira et al. 2013), and other SMGs.

Table 4. Relative Line Strengths for Average SMGs

Ratio	Value
Measurements^a
[O i]63/ $\mathrm{FIR}$	$(3.6\pm 1.2)\times {10}^{-4}$
[S iii]33/ $\mathrm{FIR}$	$\lt 3.6\times {10}^{-4}$
[Si ii]34/ $\mathrm{FIR}$	$(8.4\pm 1.7)\times {10}^{-4}$
[O iii]52/ $\mathrm{FIR}$	$\lt 3.3\times {10}^{-4}$
[N iii]57/ $\mathrm{FIR}$	$(2.7\pm 1.0)\times {10}^{-4}$
[O iv]26/ $\mathrm{FIR}$	$\lt 1.4\times {10}^{-3}$
[Fe ii]26/ $\mathrm{FIR}$	$\lt 1.4\times {10}^{-3}$
${{\rm{H}}}_{2}$ S(0)/ $\mathrm{FIR}$	$\lt 7.5\times {10}^{-4}$
${{\rm{H}}}_{2}$ S(1)/ $\mathrm{FIR}$	$\lt 1.4\times {10}^{-3}$
[C ii]158/ $\mathrm{FIR}$	$(1.7\pm 1.1)\times {10}^{-3}$
[O i]63/[C ii]158	2.2 ± 1.5

Predictions from the PDR model^b
[C i]609/ $\mathrm{FIR}$	(0.01–20) × 10⁻⁵
[C i]370/ $\mathrm{FIR}$	(0.01–60) × 10⁻⁵
[O i]145/ $\mathrm{FIR}$	(0.01–20) × 10⁻⁵
[Fe ii]26/ $\mathrm{FIR}$ , $Z={Z}_{\odot }$	(0.3–50) × 10⁻⁷
[Fe ii]26/ $\mathrm{FIR}$ , $Z=3{Z}_{\odot }$	(0.9–800) × 10⁻⁷

Notes. Upper limits are $3\sigma$ limits. $\mathrm{FIR}$ refers to the 42.5–122.5 μm continuum flux (Section 3.1).

^aFrom the mean stack measurements (Section 3.3). ^bAs discussed in Section 3.7.2, we use the best-fit parameters of the PDR model to predict the average strengths of other transitions from the PDRs of the SMGs.

Download table as: ASCII Typeset image

Table 5. Archival Sources and Observations Included in the Stacking Analyses

Name	R.A.	Decl.	z	${L}_{\mathrm{FIR}}$ ^a	FIR^a	Magnification	References^b	Program ID	OBSIDs^c
				(10¹³ ${L}_{\odot }$ )	( ${10}^{-15}\,{\mathrm{Wm}}^{-2}$ )
ID9	${09}^{{\rm{h}}}{07}^{{\rm{m}}}40\buildrel{\rm{s}}\over{.} 0$	−00°42'01''	1.577	4.4	8.2	8.8 ± 2.2	N14	OT1_averma_1	134223228[5–7], 134224524[0–2]
ID11	${09}^{{\rm{h}}}{10}^{{\rm{m}}}43\buildrel{\rm{s}}\over{.} 1$	−00°03'24''	1.786	5.7	7.7	10.9 ± 1.3	N14	OT1_averma_1	134223129[1–4]
ID17	${09}^{{\rm{h}}}{03}^{{\rm{m}}}03\buildrel{\rm{s}}\over{.} 0$	−01°41'27''	2.305	6.8	5.0	4.9 ± 0.7	N14	OT1_averma_1	134223131[3, 4]
HLock01	${10}^{{\rm{h}}}{57}^{{\rm{m}}}51\buildrel{\rm{s}}\over{.} 2$	+57°30'28''	2.958	11	4.3	10.9 ± 0.7	Co11, R11, S11	OT1_averma_1	1342232311, 134224564[4–6],
									1342256251, 1342256261
G15v2.779^d	${14}^{{\rm{h}}}{24}^{{\rm{m}}}14\buildrel{\rm{s}}\over{.} 0$	+02°23'05''	4.243	7.6	1.2	4.1 ± 0.2	B12, B13	OT1_averma_1	1342238160, 1342261470
SMM J2135^e	${21}^{{\rm{h}}}{35}^{{\rm{m}}}11\buildrel{\rm{s}}\over{.} 6$	+01°02'52''	2.326	4.3	3.6	32.5 ± 4.5	I10	OT1_averma_1	1342231704, 1342244443, 1342245235,
									1342245393, 134225694[0–1], 1342257256
NC.v1.143	${12}^{{\rm{h}}}{56}^{{\rm{m}}}32\buildrel{\rm{s}}\over{.} 6$	+23°36'27''	3.565	8.9	2.2	11.3 ± 1.7	B13, Rp	OT1_averma_1	1342257[799, 800]
NA.v1.177	${13}^{{\rm{h}}}{28}^{{\rm{m}}}59\buildrel{\rm{s}}\over{.} 3$	+29°23'27''	2.778	6.0	2.7	...	B13	OT1_averma_1	134225956[0, 1]
SWIRE 3-9	${10}^{{\rm{h}}}{43}^{{\rm{m}}}43\buildrel{\rm{s}}\over{.} 9$	+57°13'23''	1.735	0.47	0.73	1	B15	OT2_dbrisbin_1	1342253586, 1342253776
SWIRE 3-14	${10}^{{\rm{h}}}{45}^{{\rm{m}}}14\buildrel{\rm{s}}\over{.} 5$	+57°57'09''	1.780	0.28	0.29	1	B15	OT2_dbrisbin_1	1342247014, 1342247131
SWIRE 4-5	${10}^{{\rm{h}}}{44}^{{\rm{m}}}27\buildrel{\rm{s}}\over{.} 5$	+58°43'10''	1.756	0.10	0.09	1	B15	OT2_dbrisbin_1	134224663[8, 9]
SWIRE 4-15	${10}^{{\rm{h}}}{46}^{{\rm{m}}}56\buildrel{\rm{s}}\over{.} 5$	+59°02'36''	1.854	0.30	0.37	1	B15	OT2_dbrisbin_1	1342253587, 1342253775
SDSS J1206	${12}^{{\rm{h}}}{06}^{{\rm{m}}}01\buildrel{\rm{s}}\over{.} 7$	+51°42'28''	1.999	0.42	0.46	∼27	B15	OT2_dbrisbin_1	1342246801
SMM J0302	${03}^{{\rm{h}}}{02}^{{\rm{m}}}27\buildrel{\rm{s}}\over{.} 7$	+00°06'52''	1.408	0.46	1.2	1	B15	OT2_dbrisbin_1	134224778[4, 5]
MIPS 22530	${17}^{{\rm{h}}}{23}^{{\rm{m}}}03\buildrel{\rm{s}}\over{.} 3$	+59°16'00''	1.950	0.62	0.69	1	B15	OT2_dbrisbin_1	1342249495, 1342256260
LESS21	${03}^{{\rm{h}}}{33}^{{\rm{m}}}29\buildrel{\rm{s}}\over{.} 7$	−27°34'44''	1.235	0.07	0.13	1	C12	OT1_kcoppin_1	1342239701
LESS34	${03}^{{\rm{h}}}{32}^{{\rm{m}}}17\buildrel{\rm{s}}\over{.} 6$	−27°52'28''	1.098	0.06	0.14	1	C12	OT1_kcoppin_1	1342239703
LESS66	${03}^{{\rm{h}}}{33}^{{\rm{m}}}31\buildrel{\rm{s}}\over{.} 9$	−27°54'10''	1.315	0.14	0.41	1	C12	OT1_kcoppin_1	1342239369
LESS88	${03}^{{\rm{h}}}{31}^{{\rm{m}}}54\buildrel{\rm{s}}\over{.} 8$	−27°53'41''	1.269	0.08	0.17	1	C12	OT1_kcoppin_1	1342239705
LESS106	${03}^{{\rm{h}}}{31}^{{\rm{m}}}40\buildrel{\rm{s}}\over{.} 2$	−27°56'22''	1.617	0.15	0.23	1	C12	OT1_kcoppin_1	1342239753
LESS114	${03}^{{\rm{h}}}{31}^{{\rm{m}}}51\buildrel{\rm{s}}\over{.} 1$	−27°44'36''	1.606	0.33	0.57	1	C12	OT1_kcoppin_1	1342239702
SPT 0538-50	${05}^{{\rm{h}}}{38}^{{\rm{m}}}16\buildrel{\rm{s}}\over{.} 5$	−50°30'50''	2.782	5.1	2.3	21.0 ± 4.0	Bo13	OT2_dmarrone_2	1342270691
SPT 0125-47	${01}^{{\rm{h}}}{25}^{{\rm{m}}}07\buildrel{\rm{s}}\over{.} 0$	−47°23'57''	2.515	8.7	4.9	5.5 ± 0.1	A16, V13, We13	OT2_dmarrone_2	1342270768
SPT 0103-45	${01}^{{\rm{h}}}{03}^{{\rm{m}}}11\buildrel{\rm{s}}\over{.} 4$	−45°38'54''	3.092	3.4	1.2	7.2 ± 5.2	G15, V13, We13, S16	OT2_dmarrone_2	1342271050
F10214	${10}^{{\rm{h}}}{24}^{{\rm{m}}}34\buildrel{\rm{s}}\over{.} 6$	+47°09'09''	2.286	8.4	6.0	∼12	S10	SDP_kmeisenh_3	1342186812, 1342187021
SMM J14011	${14}^{{\rm{h}}}{01}^{{\rm{m}}}04\buildrel{\rm{s}}\over{.} 9$	+02°52'24''	2.565	1.4	0.76	3.5 ± 0.5	S05, S13	KPGT_kmeisenh_1	134221331[1–4], 1342213677
SMM J22471	${22}^{{\rm{h}}}{47}^{{\rm{m}}}10\buildrel{\rm{s}}\over{.} 4$	−02°05'53''	1.158	2.7	10	∼2	S10	OT1_gstacey_3	1342211842, 1342212211
SWIRE J104738	${10}^{{\rm{h}}}{47}^{{\rm{m}}}38\buildrel{\rm{s}}\over{.} 3$	+59°10'10''	1.958	0.40	0.42	1	S10	OT1_gstacey_3	134223226[8, 9]
SWIRE J104704	${10}^{{\rm{h}}}{47}^{{\rm{m}}}05\buildrel{\rm{s}}\over{.} 1$	+59°23'33''	1.954	1.0	1.1	1	S10	OT1_gstacey_3	134223227[0, 1]
SMM J123634	${12}^{{\rm{h}}}{36}^{{\rm{m}}}34\buildrel{\rm{s}}\over{.} 6$	+62°12'41''	1.222	0.54	1.8	1	S10	OT1_gstacey_3	1342232[599, 601]
MIPS J142824	${14}^{{\rm{h}}}{28}^{{\rm{m}}}24\buildrel{\rm{s}}\over{.} 1$	+35°26'18''	1.325	1.3	3.6	$\lt 10$	HD10, S10	SDP_esturm_3	1342187779
SMM J02396	${02}^{{\rm{h}}}{39}^{{\rm{m}}}56\buildrel{\rm{s}}\over{.} 7$	−01°34'24''	1.062	0.17	0.84	∼2.3	G05, C11, C13	KPGT_esturm_1	1342214674

Notes.

^a ${L}_{\mathrm{FIR}}$ ${L}_{\mathrm{FIR}}$ and FIR are the apparent (i.e., not corrected for lensing) far-IR luminosity (40–500 μm) and continuum flux (42.5–122.5 μm), respectively (Section 3.1). Where unavailable from published SED fits, these are estimated by scaling published infrared luminosities to the required rest-frame wavelength ranges using the average fitted SED of the targeted SMGs. ^bReferences are not the complete list of published papers on each source, but rather the source of data used in this paper: A16: Aravena et al. (2016), Bo13: Bothwell et al. (2013), B12: Bussmann et al. (2012), B13: Bussmann et al. (2013), B15: Brisbin et al. (2015), C11: Chen et al. (2011), Co11: Conley et al. (2011), C12: Coppin et al. (2012), C13: Chen et al. (2013), G05: Greve et al. (2005), G15: Gullberg et al. (2015), HD10: Hailey-Dunsheath et al. (2010), I10: Ivison et al. (2010b), N14: Negrello et al. (2014), R11: Riechers et al. (2011b), Rp: D. Riechers et al. (2017, in preparation), S05: Smail et al. (2005), S10: Stacey et al. (2010), S11: Scott et al. (2011), S13: Sharon et al. (2013), S16: Spilker et al. (2016), V13: Vieira et al. (2013), We13: Weis et al. (2013). ^cOBSIDs are the Herschel observation identification number(s) for this program, used to identify the the photometric and spectroscopic observation of each target in the Herschel archive. ^dAlso known as ID15.141 or ID141. ^eAlso known as the Cosmic Eyelash.

Download table as: ASCII Typeset image

The PACS spectroscopy of the archival targets is reduced in the same way as the targeted data (Section 2.2). For those spectra that have been published elsewhere, we have verified that our reduction produces measurements consistent with the published data. PACS photometry is not available for most of the archival targets, so those are not considered here; we instead use the published IR luminosity of each source. Where necessary, we scale to the wavelength ranges for ${L}_{\mathrm{FIR}}$ and ${FIR}$ (Section 3.1) by using the SED fits of the targeted sources (Section 3.1). For sources with multiple published IR luminosities, we use the one constrained by the most photometric data points.

3. Analysis and Discussion

3.1. Far-Infrared SED Fits

The PACS photometry, described in Section 2.3, is supplemented with the SPIRE (Griffin et al. 2010) 250, 350, and 500 μm data from HerMES (Roseboom et al. 2012; Wang et al. 2014) and H-ATLAS (Valiante et al. 2016), and, where available, longer wavelength follow-up photometry (see references in Table 1). We show the far-IR SEDs derived from this compilation of data in Figure 1.

**Figure 1.** Far-IR SEDs of the 13 targets in our sample, combining existing photometry (see references in Table 1) with our new *Herschel*-PACS 70 and 160 μm data and 100 μm measurements from R. George et al. (2017, in preparation). The best-fit modified black body SEDs, which are used to calculate ${FIR}$ and $\mathrm{far-IR\; continuum\; flux}$ (see Section 3.1) are shown. For each galaxy, we also present the best-fit SED from the Dale & Helou (2002) template library.
Download figure:
Standard image High-resolution image

**Figure 1.** Far-IR SEDs of the 13 targets in our sample, combining existing photometry (see references in Table 1) with our new *Herschel*-PACS 70 and 160 μm data and 100 μm measurements from R. George et al. (2017, in preparation). The best-fit modified black body SEDs, which are used to calculate ${FIR}$ and $\mathrm{far-IR\; continuum\; flux}$ (see Section 3.1) are shown. For each galaxy, we also present the best-fit SED from the Dale & Helou (2002) template library.
Download figure:
Standard image High-resolution image

For each galaxy, we fit the observed far-infrared SED with an optically thin modified blackbody spectrum of the form

$\begin{eqnarray}&&{S}_{\nu }\propto {\nu }^{\beta }{B}_{\nu }({T}_{{\rm{D}}}),\end{eqnarray} \tag{ 1 }$

where ${S}_{\nu }$ is the flux density, ν is frequency, β is the power law emissivity index, and ${B}_{\nu }({T}_{{\rm{D}}})$ is the Planck function, defined as

$\begin{eqnarray}&&{B}_{\nu }({T}_{{\rm{D}}})\propto \displaystyle \frac{{\nu }^{3}}{{e}^{h\nu /{{kT}}_{{\rm{D}}}}-1},\end{eqnarray} \tag{ 2 }$

for a dust temperature, ${T}_{{\rm{D}}}$ . Here, h and k denote the Planck and Boltzmann constants, respectively. We fix $\beta =1.5$ , which is consistent with observed values in a range of galaxies (e.g., Hildebrand 1983; Dunne & Eales 2001), and allow ${T}_{{\rm{D}}}$ and the normalization to vary. The best-fit modified blackbody curve for each galaxy is shown in Figure 1.

Using these modified blackbody fits, we next calculate both far-IR luminosity ( ${L}_{\mathrm{FIR}}$ ) and $\mathrm{far-IR\; continuum\; flux}$ ( ${FIR}$ ) for each SMG. For consistency with existing studies, we follow the definitions of Graciá-Carpio et al. (2011) and Coppin et al. (2012) for these quantities, whereby:

1.
${L}_{\mathrm{FIR}}$ is the luminosity of the rest-frame SED integrated between 40 and 500 μm, and;
2.
${FIR}$ is the luminosity integrated between 42.5 and 122.5 μm in the rest-frame, and converted to flux by dividing by $4\pi {D}_{{\rm{L}}}^{2}$ , where ${D}_{{\rm{L}}}$ is the luminosity distance.

The apparent (i.e., without correction for lensing amplification) values of ${L}_{\mathrm{FIR}}$ and ${FIR}$ , calculated from the modified blackbody SED fits, are listed in Section 3.1 and are used in the analysis in the rest of this paper.

However, because the single temperature modified blackbody can underpredict the emission on the Wien side of the far-IR dust peak, it is possible that the ${L}_{\mathrm{FIR}}$ and ${FIR}$ values that we calculate from the modified blackbody fits are systematically underestimated. To test the magnitude of this effect, we also fit each galaxy with SEDs from the Dale & Helou (2002) template library; these fits are also shown in Figure 1.

There is no significant systematic offset between ${L}_{\mathrm{FIR}}$ from the two fitting methods, with the median ratio of the Dale & Helou (2002) to modified blackbody values being 0.99. There are only three galaxies with ${L}_{\mathrm{FIR}}$ from the Dale & Helou (2002) SED fits that differ significantly from the values of the modified blackbody fits. These are G15v2.19, G12v2.43, and HBoötes01, the first of which is 15% higher, and both of the latter are 10% lower for the Dale & Helou (2002) fits. Only one galaxy has significantly higher ${FIR}$ from the Dale & Helou (2002) SEDs than the modified blackbody fits: G12v2.257, a with 30% difference in ${FIR}$ . However, five systems—G15v2.19, G09v1.40, NGP.NA.144, HXMM01, and NGP.NB.78—have lower ${FIR}$ for the Dale & Helou (2002) SEDs than the modified blackbody fits. For these six galaxies, the ${FIR}$ for the Dale & Helou (2002) fits are 70%–95% of the modified blackbody values. These galaxies would thus be offset upward by $\sim 25 \%$ in Figures 3 and 4 if we were to use ${FIR}$ from the Dale & Helou (2002) SED fits instead of from the modified blackbody fits. The typically small differences between ${L}_{\mathrm{FIR}}$ and ${FIR}$ from the Dale & Helou (2002) and modified blackbody fits are because ${L}_{\mathrm{FIR}}$ and ${FIR}$ are most sensitive to the peak and long wavelength part of the SED, where the modified blackbody does a good job of fitting the data. Note that, due to the narrow wavelength ranges considered for ${FIR}$ and ${L}_{\mathrm{FIR}}$ , and the slightly lower normalization of some of the Dale & Helou (2002) fits, having lower ${L}_{\mathrm{FIR}}$ for the Dale & Helou (2002) fits compared with the modified blackbody in some cases is not unexpected.

**Figure 2.** Mean stacks of the continuum-subtracted rest-frame spectra for each of the observed IR emission lines, with the best-fit $\geqslant 3\sigma$ Gaussian line profiles overlaid (red). The lower (blue) line corresponds to the right-hand axis and shows the number of spectra included in each bin of the stack, which is variable because of the different rest-frame coverage of each observation.
Download figure:
Standard image High-resolution image

$\geqslant 3\sigma $ — **Figure 2.** Mean stacks of the continuum-subtracted rest-frame spectra for each of the observed IR emission lines, with the best-fit $\geqslant 3\sigma$ Gaussian line profiles overlaid (red). The lower (blue) line corresponds to the right-hand axis and shows the number of spectra included in each bin of the stack, which is variable because of the different rest-frame coverage of each observation.
Download figure:
Standard image High-resolution image

**Figure 3.** Emission line to continuum ${FIR}$ flux (42.5–122.5 μm) as a function of (lensing corrected) ${L}_{\mathrm{FIR}}$ (40–500 μm) for the observed fine structure lines. Our individual SMG targets are shown in black, with $3\sigma$ upper limits plotted for transitions that are not measured above this level. For individual observations, error bars represent $1\sigma$ uncertainties, including contributions from the line fitting, from the ${L}_{\mathrm{FIR}}$ and ${FIR}$ measurements, and from the lensing amplification. The star represents results from mean stacking (Section 3.3), using the average ${L}_{\mathrm{FIR}}$ and ${FIR}$ of the galaxies included in each stack. For comparison, colored symbols show a compilation of star-forming galaxies, AGN, LINERs, and unclassified galaxies (see Section 3.4), with the filled symbols representing those at $z\gt 0.5$ (Farrah et al. 2007; O'Halloran et al. 2008; Bernard-Salas et al. 2009; Sturm et al. 2010; Graciá-Carpio et al. 2011; Coppin et al. 2012).
Download figure:
Standard image High-resolution image

**Figure 3.** Emission line to continuum ${FIR}$ flux (42.5–122.5 μm) as a function of (lensing corrected) ${L}_{\mathrm{FIR}}$ (40–500 μm) for the observed fine structure lines. Our individual SMG targets are shown in black, with $3\sigma$ upper limits plotted for transitions that are not measured above this level. For individual observations, error bars represent $1\sigma$ uncertainties, including contributions from the line fitting, from the ${L}_{\mathrm{FIR}}$ and ${FIR}$ measurements, and from the lensing amplification. The star represents results from mean stacking (Section 3.3), using the average ${L}_{\mathrm{FIR}}$ and ${FIR}$ of the galaxies included in each stack. For comparison, colored symbols show a compilation of star-forming galaxies, AGN, LINERs, and unclassified galaxies (see Section 3.4), with the filled symbols representing those at $z\gt 0.5$ (Farrah et al. 2007; O'Halloran et al. 2008; Bernard-Salas et al. 2009; Sturm et al. 2010; Graciá-Carpio et al. 2011; Coppin et al. 2012).
Download figure:
Standard image High-resolution image

**Figure 4.** Same as Figure 3, for the H₂ rotational transitions.
Download figure:
Standard image High-resolution image

3.2. Individual Emission Line Measurements

We measure the line fluxes (and upper limits) from the 1D spectra of the individual galaxies, reduced and extracted as described in Section 2.2, with the noisy regions at the edges of the spectra (typically 5–10 wavelength bins) removed. Then, with the exception of the [O iv]26 observations, each spectrum is fit with a single Gaussian line profile and flat continuum component using the mpfit function in IDL (Markwardt 2009), which uses nonlinear Levenberg-Marquardt least-squares minimization. We constrain the fits to have non-negative continua. The velocity offsets of the lines are required to be $\leqslant \pm 800\,{\mathrm{kms}}^{-1}$ from their expected locations, based on the CO redshifts. The wavelength range of the [O iv]26 observations includes the [Fe ii]26 line. Therefore, those data are fitted with double Gaussians, using the same mpfit IDL function. In all cases, the velocity-integrated flux in each line is calculated from the continuum-subtracted best-fit Gaussian.

The pipeline-derived uncertainties on the PACS spectra are known to be unreliable.²⁵ Therefore, we weigh each wavelength bin equally for fitting purposes. The uncertainty on the line fluxes are determined from 1000 trials for each line, wherein we add random noise with the same $1\sigma$ rms as measured from the line-free portions of the spectra and refit the line. The $3\sigma$ detection limit for each line is calculated from a Gaussian profile with a peak height three times the rms noise in the spectra, centered at the expected position of the emission line from the CO redshift. For the purposes of this calculation, we assume a linewidth of 300 km s⁻¹ FWHM, which is consistent with observations of high-redshift star-forming galaxies (e.g., Sturm et al. 2010; Coppin et al. 2012) and similar to the PACS instrumental resolution.

The spectra and line fits for the 13 targets of OT2_jwardlow_1 are presented in Appendix B, and the measurements given in Table 3.

3.3. Stacked Spectra

We next investigate the average properties of the spectra by stacking the observations of each transition for all the galaxies, which reduces the background noise by a factor of $\sim \sqrt{N}$ for a stack of N galaxies with the same background. To trace to fainter noise limits, we include both the 13 targeted galaxies, and the 32 archival sources in the stacks. We have verified that the measured line fluxes (or limits) are consistent, whether the archival data are included or not. For each line, the stacked spectra contains 8–37 galaxies. Therefore, we expect improvements of factors of ∼3–6 in the average sensitivity of individual spectra by stacking.

To perform the stacking, we first shift each spectrum to the rest frame and subtract the continuum. We generate a base rest-frame wavelength grid with spacing equal to the average rest-frame native PACS resolution for each line targeted. The individual spectra are then rebinned to the new rest-frame wavelength grid and three different stacks are generated.

1.
Our fiducial method is a mean stack, derived by calculating the mean value in each wavelength bin. These mean spectra for each targeted emission line are shown in Figure 2, and the measurements are presented in Table 3. For the line flux/ ${FIR}$ ratios examined in Sections 3.4 and 3.7, we use the mean ${FIR}$ (42.5–122.5 μm) of the sources included in the stack, such that the ratio is equivalent to mean(line flux)/mean( ${FIR}$ ).
2.
We also generate median stacks, consisting of the median value in each wavelength bin, which are used to investigate whether a few bright outliers dominate the fiducial mean stacks.
3.
To investigate the presence of trends with infrared emission, weighted mean stacks are also produced, where each source is weighted by $1/{FIR}$ (42.5–122.5 μm), In this case, measurements from the weighted mean stacks are equivalent to mean(line flux/ ${FIR}$ ).

The rest-frame wavelength coverage from different observing programs varies, so the number of galaxies contributing to each wavelength bin varies, as is shown in Figure 2. In each case, the stacked spectra are fit using the methodology described in Section 3.2 for the individual observations. Because the number of data points stacked in each wavelength bin varies, the noise level is weighted across the spectra according to $1/\sqrt{N}$ . This is valid because the observations for each line have similar depths. For the fiducial mean stacks, the measured average fluxes (or limits) are reported in Table 3, and the line/ ${FIR}$ (42.5–122.5 μm) ratios in Table 4.

The measured linewidths for the [O i]63, [Si ii]34, and [N iii]57 (those with $\geqslant 3\sigma$ detections) in the fiducial mean stacks are 510 ± 160, 820 ± 280, and $700\pm 240\,{\mathrm{kms}}^{-1}$ , respectively (Figure 2). We caution that these apparent linewidths are not physically meaningful, because they include a contribution from for potential offsets between the literature spectroscopic redshifts (due to broadband CO searches or optical data; Sections 2.1 and 2.4) and the targeted IR transitions, which will artificially broaden the lines.

The line fluxes and upper limits are consistent between the median and (fiducial) mean stacks, which demonstrates that the mean stacks are not dominated by a few bright outliers. With the exception of the [O i]63 line, the weighted stack measurements are also consistent with the mean stacks, showing that, for most of the lines, there is no evidence of correlations with infrared emission for SMGs. For [O i]63, there is no detection in the weighted stack, with a $3\sigma$ detection limit of [O i]/ ${FIR}$ $\lt \,3\times {10}^{-3}$ (compared with [O i]/ ${FIR}$ = $(3.6\pm 1.2)\times {10}^{-4}$ for the fiducial mean stack). This suggests that, for SMGs, there may be an inverse correlation between infrared luminosity and [O i]63 emission. The rest of this paper focuses on the fiducial mean stacked fluxes, but we discuss, where relevant, how the conclusions would change if we instead considered the [O i]63 weighted stack measurement.

3.4. Individual Line Strengths

One way to characterize the strength of IR emission lines is via the line to ${FIR}$ (42.5–122.5 μm) ratio (e.g., Fischer et al. 2010; Sturm et al. 2010; Graciá-Carpio et al. 2011; Coppin et al. 2012; Magdis et al. 2014). These are presented in Figures 3 and 4 for the fine structure lines and ${{\rm{H}}}_{2}$ lines, respectively, and discussed here. Measurements of the mean stacked spectroscopy for our lensed SMGs are presented, with the relevant ${L}_{\mathrm{FIR}}$ and ${FIR}$ calculated, as the mean of the galaxies included in the stack. These derived average line flux-to- $\mathrm{FIR}$ ratios for SMGs are given in Table 4. Published measurements for other galaxies (mostly at low redshift) are also shown in Figures 3 and 4, color-coded by whether they are star-forming galaxies, AGN, LINERs, or unclassified (Colbert et al. 1999; Malhotra et al. 2001; Negishi et al. 2001; Sturm et al. 2002, 2006, 2010; Lutz et al. 2003; Verma et al. 2003; Dale et al. 2004; Farrah et al. 2007, 2013; Brauher et al. 2008; O'Halloran et al. 2008; Tommasin et al. 2008, 2010; Bernard-Salas et al. 2009; Hao et al. 2009; Veilleux et al. 2009; Hunt et al. 2010; Ivison et al. 2010a; Graciá-Carpio et al. 2011; Valtchanov et al. 2011; Coppin et al. 2012; Stierwalt et al. 2014). Note that most of the targeted emission lines ([O i]63, [S iii]33, [O iii]52 and [N iii]57) predominantly trace PDRs and H ii regions. Any weak AGN contribution will decrease the relative line to ${FIR}$ ratio, as the continuum emission is preferentially enhanced. Although energetically dominant or very powerful AGN can sometimes contribute to the line flux, such AGN are exceptionally rare in SMGs (e.g., Sections 3.4.6 and 3.5; Alexander et al. 2005; Valiante et al. 2007; Pope et al. 2008; Menéndez-Delmestre et al. 2009; Laird et al. 2010; Wang et al. 2013) and are unlikely to affect our measurements.

Locally, the relative strength of many PDR cooling lines, including [O i]63, [S iii]33, [N ii]122, and [C ii]158, are suppressed with respect to the far-IR emission in the most luminous systems, particularly those with "warmer" infrared colors (e.g., Malhotra et al. 2001; Graciá-Carpio et al. 2011; Farrah et al. 2013). Various explanations for the emission line deficits in high-luminosity galaxies have been proposed, including their being dustier and having higher ionization parameters in the ISM, resulting in a higher fraction of the UV photons being absorbed by dust and re-emitted in the far-IR, enhancing the far-IR brightness, and thus decreasing the line/FIR ratios (e.g., Luhman et al. 2003; González-Alfonso et al. 2008; Abel et al. 2009; Graciá-Carpio et al. 2011; Farrah et al. 2013; Fischer et al. 2014). Alternative explanations include non-PDR flux in the far-IR, such as from AGN, which would also serve to dilute the PDR line emission (e.g., Malhotra et al. 2001; Luhman et al. 2003; Farrah et al. 2013), or the primary gas coolant not being the typical [C ii]158 or [O i]63 lines, but instead via other mechanisms (e.g., Farrah et al. 2013). Therefore, it is probable that PDR line deficits may be indicative of a different "mode" of star formation in local ULIRGs compared with sub-LIRGs, with the ULIRGs' star formation being more concentrated, as is typical in merger-induced activity. We next probe whether SMGs exhibit similar deficits on a transition-by-transition basis.

3.4.1. [O i]63 $\mu {\rm{m}}$

[O i]63 has a critical density of $\sim 5\times {10}^{5}\,{\mathrm{cm}}^{-3}$ , traces dense molecular gas with ${\rm{T}}\gt 100$ K, and is one of the dominant cooling lines in dense PDRs. As can be seen in Figure 3, local LIRGs and ULIRGs have an [O i]63 deficit compared with lower luminosity systems, which typically have [O i]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}\sim 2\,\times \,{10}^{-3}$ (e.g., Luhman et al. 2003; Graciá-Carpio et al. 2011; Farrah et al. 2013).

Initial observations and publications hinted that [O i]/ $\mathrm{FIR}$ may be enhanced in SMGs, with [O i]/ $\mathrm{FIR}$ similar to local sub-LIRGs (Sturm et al. 2010; Coppin et al. 2012), implying that they may have large reservoirs of dense gas (to fuel the far-IR luminosities, and probed with CO) like local ULIRGs, but with star formation efficiencies comparable to late-type galaxies. However, the additional data from our observations (Figure 3) now show that the picture is more complicated, with only four out of all 15 SMGs ever observed in [O i]63 (MIPS J142824.0+352619; Sturm et al. 2010, and unlensed examples from Coppin et al. 2012) detected (although mostly at low significance). The flux limits available for the majority of the remainder of individual galaxies are not deep enough to provide robust constraints, leaving interpretation of those results open to discussion.

Our mean stacked data are significantly more constraining, providing a $3.2\sigma$ detection, with [O i]/ $\mathrm{FIR}$ = $(0.36\pm 0.12)\times {10}^{-3}$ . Thus, the mean stack result indicates that, on average, high-redshift SMGs behave like local ULIRGs, with a deficit in their [O i]63 emission relative to ${FIR}$ , although there are exceptions. The non-detection of [O i]63 in the ${FIR}$ -weighted stack also suggests that there may be a trend in the strength of [O i]63 emission with ${FIR}$ for SMGs. The difference between the [O i]/ $\mathrm{FIR}$ in the few individually detected SMGs (from Sturm et al. 2010 and Coppin et al. 2012) and the mean stack is significant, suggesting there may be physical differences between them, perhaps including a range of possible [O i]63 and ${FIR}$ emission mechanisms for the SMGs. The two galaxies with the most compelling detections in Coppin et al. (2012) both potentially contain AGN (see their discussion), which may be an explanation, as weak AGN can strengthen the relative line flux (notice that the local sources with the highest [O i]/ $\mathrm{FIR}$ are AGN and LINERs). We further investigate the [O i]63 emission in Section 3.7, where we use PDR models to probe the state of the ISM in high-redshift SMGs.

3.4.2. [S iii]33 $\mu {\rm{m}}$

[S iii]33 has a lower critical density than [O i]63, and is a key coolant of H ii regions. Figure 3 shows that the [S iii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ ratio is locally anti-correlated with ${L}_{\mathrm{FIR}}$ (40–500 μm) for star-forming galaxies and AGN. This may be an effect of continuum dilution of the line, as is observed in other fine-structure lines (e.g., [C ii]158; Dale et al. 2006). Neither our stack nor any of our individual SMGs have $\geqslant 3\sigma$ detections in [S iii]33. The upper limits for the undetected individual SMGs are consistent with luminous local galaxies. The limit on the average [S iii]/FIR of SMGs from the stacked data is lower than many local systems, but consistent with expectations if the local trend is extrapolated to the higher luminosity of the stack. Our stacked data are also consistent with the one other observation of [S iii]33 in a high-redshift source (IRASF10214+4724; Sturm et al. 2010).

3.4.3. [Si ii]34 $\mu {\rm{m}}$

[Si ii]34 is an important cooling line in PDRs. It has a higher ionization energy than [O i]63, but a similar critical density; thus, higher-intensity radiation fields are required to excite [Si ii]34 compared to [O i]63, and it can also be emitted from XDRs. Observations of [Si ii]34 have so far been limited to the local universe, where [Si ii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}\sim {10}^{-3}$ for AGN, star-forming galaxies, and LINERs. Locally, there is hint of a trend of lower [Si ii]/ $\mathrm{FIR}$ in the highest luminosity systems, but not as convincingly as for [O i]63 or [C ii]158 (for example). There are also indications that the highest [Si ii]/ $\mathrm{FIR}$ values (≳4 × 10⁻³) are only present in LINERs and AGN, i.e., high ionization environments (Section 3.5). None of the 11 individual SMGs targeted here are detected at the $\geqslant 3\sigma$ level. However, the [Si ii]34 stack is our strongest detection of all the data (7.5σ), and shows that, on average, SMGs have [Si ii]/ $\mathrm{FIR}$ = $(8.4\pm 1.7)\times {10}^{-4}$ . This average value of [Si ii]/ $\mathrm{FIR}$ for SMGs is marginally lower than the local average for all galaxies, but substantially higher than may be expected if the weak trend of [Si ii]34/ $\mathrm{FIR}$ dropping for high-luminosity local galaxies continues into the ULIRG regime.

3.4.4. [O iii]52 $\mu {\rm{m}}$

[O iii]52 is an efficient tracer of H ii regions, and local ULIRGs show the same suppression in [O iii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ as for many of the other fine-structure lines. In our [O iii]52 observations, one SMG is detected at $\gt 3\sigma$ (NGP.NA.144); the remaining seven individually targeted SMGs and the stack are all undetected. In all cases, our constraints on [O iii]/ $\mathrm{FIR}$ are consistent with the local LIRG and sub-LIRG population. However, in most cases, we cannot rule out lower ratios for SMGs (e.g., as is observed in F10214 at $z=2.28;$ Sturm et al. 2010). The detected SMG, NGP.NA.144, has [O iii]/ $\mathrm{FIR}$ = $(3.4\pm 1.0)\times {10}^{-3}$ , which is significantly higher then local ULIRGs, but is somewhat consistent with local LIRGs and sub-LIRGs. In contrast, the stack has [O iii]/ $\mathrm{FIR}$ $\lt \,0.33\times {10}^{-3}$ , consistent with local ULIRGs and lower than most detected local star-forming galaxies. Similarly to the [O i]63 data (Section 3.4.1), the dichotomy between the individual NGP.NA.144 detection and the stacked [O iii]52 measurement is suggestive of a range of conditions in the [O iii]52 emitting region of SMGs.

3.4.5. [N iii]57 $\mu {\rm{m}}$

[N iii]57 has similar ionization energy and critical density to [O iii]52, and is also a strong tracer of H ii regions. Local ULIRGs have a similar [N iii]57 deficit to [O iii]52, and other fine structure lines. We targeted [N iii]57 in four lensed SMGs, and none were individually detected above $3\sigma$ , although the [N iii]57 line is detected in the stack of eight SMGs at $3.0\sigma$ . The stack, as well as the four individual targets, have [N iii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ values that are broadly consistent with local sub-LIRGs, although still within the upper range of values of local ULIRGs. Overall, the [N iii]57 emission regions in SMGs likely have similar conditions to local star-forming galaxies and some local ULIRGs.

3.4.6. [O iv]26 $\mu {\rm{m}}$

[O iv]26 is a high-excitation line. As such, it is a reliable tracer of AGN activity (e.g., Meléndez et al. 2008; Rigby et al. 2009). As can be seen in Figure 3, all local sources with [O iv]/FIR ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}\gtrsim {10}^{-3}$ are AGN. Other galaxy types can have some [O iv]26 emission, but it is always fainter (relative to ${FIR}$ ) than this limit, because the [O iv]26 contribution is enhanced in AGN environments.

Of our 10 [O iv]26 targets, nine are undetected and one is detected at $3.1\sigma$ (HXMM01). All of the undetected systems have upper limits of [O iv]/ $\mathrm{FIR}$ ≃ 2 × 10⁻². Therefore, we cannot rule out some AGN contribution to these galaxies. For HXMM01, we measure [O iv]/ $\mathrm{FIR}$ = $(3.9\pm 1.3)\times {10}^{-3}$ , which is significantly above the observed ratios in local star-forming galaxies, and an indication of a hidden AGN in this system. HXMM01 was previously studied by Fu et al. (2013), who found no evidence of an AGN in the mid-IR (IRAC) or far-IR colors, or in (shallow) X-ray observations. They did observe a broad Hα line in HXMM01, but concluded that it was likely driven by starburst outflows. Thus, our [O iv]26 observation is currently the only distinct evidence of an AGN in HXMM01.

[O iv]26 is not detected in the stacked spectrum, placing an upper limit on [O iv]/ $\mathrm{FIR}$ of $1.4\times {10}^{-3}$ in average SMGs—within the region inhabited by local AGN, star-forming galaxies, and LINERs. Thus, there is no evidence from the [O iv]26 that typical SMGs contain strong AGN, but the current data do not enable us to rule out lower-luminosity AGN. This is consistent with previous analyses of SMGs using other AGN tracers (e.g., Alexander et al. 2005; Pope et al. 2008; Wang et al. 2013), and broadly consistent with the results from [Si ii]34 and [S iii]33 (Section 3.5). This hints at some AGN emission if the AGN in SMGs are weak, and therefore not picked up by our relatively shallow [O iv]26 data.

3.4.7. [Fe ii]26 $\mu {\rm{m}}$

[Fe ii]26 was not specifically targeted by this program, as it is rarely detectable, even locally. However, it occurs only ∼0.1 μm (rest-frame) away from [O iv]26 and is therefore included in the spectral coverage of the 10 SMGs for which we targeted the [O iv]26 transition. [Fe ii]26 is predominantly emitted from PDRs, and it has a low ionization energy and high critical density, similar to [Si ii]34. Also similarly to [Si ii]34, [Fe ii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ is weakly anti-correlated with ${L}_{\mathrm{FIR}}$ (40–500 μm) locally (Figure 3), although [Fe ii]26 has been measured in fewer galaxies for comparison (e.g., Verma et al. 2003; Farrah et al. 2007; O'Halloran et al. 2008; Bernard-Salas et al. 2009). As our [Fe ii]26 data were also obtained serendipitously, they are shallower than required to probe the local observed range of [Fe ii]/ $\mathrm{FIR}$ . All the individual targets, as well as the stack, are undetected, with limits consistent with local observations.

3.4.8. H₂ S(0) and H₂ S(1)

The ${{\rm{H}}}_{2}$ S(0) and ${{\rm{H}}}_{2}$ S(1) rotational transitions trace warm ( $T\simeq 150$ K) gas reservoirs (e.g., Rigopoulou et al. 2002; Roussel et al. 2007; Nesvadba et al. 2010; Higdon et al. 2014), winds (e.g., Beirão et al. 2015) and likely shocked gas in extreme systems with enhanced ${{\rm{H}}}_{2}$ emission (e.g., Appleton et al. 2006). Thus, with detections of ${{\rm{H}}}_{2}$ S(0) and ${{\rm{H}}}_{2}$ S(1), we would be able to directly measure the warm ${{\rm{H}}}_{2}$ gas mass (rather than purely relying on CO), and use the prevalence of shocks to identify violent merger activity. However, the data are not as deep as expected (Figure 4), and neither the stack nor any of the individual SMGs have emission detectable at the $\geqslant 3\sigma$ level.

3.5. AGN Contribution Traced by [Si ii]34 $\mu {\rm{m}}$ and [S iii]33 $\mu {\rm{m}}$ Emission

Previous studies have shown that the ratio of [Si ii]34 to [S iii]33 flux is an effective discriminator between AGN, LINERs, and H ii regions (Dale et al. 2006, 2009). Typically, the [Si ii]34/[S iii]33 ratio is used in conjunction with a second discriminator such as the ratios of [Ne iii] 16 μm to [Ne ii] 13 μm (Dale et al. 2006) or [Fe ii] 26 μm to [Ne ii] 13 μm (Dale et al. 2009). Observations of these additional lines are not available for SMGs, but the [Si ii]34/[S iii]33 ratio alone is still useful.

Five of our individual targets and the stack have data for both the [Si ii]34 and [S iii]33 transitions. However, none have $\geqslant 3\sigma$ detections in one or both of these lines. Therefore, the AGN contribution cannot be traced on a individual target basis from our [Si ii]34 and [S iii]33 data. For the stacked data, the [Si ii]34 detection and the [S iii]33 limit give [Si ii]34/[S iii]33 $\geqslant 3.35$ ( $3\sigma$ ), placing the average SMG in region I and II of Dale et al. (2006), corresponding to AGN and LINER emission. Thus, the [Si ii]34 and [S iii]33 measurements suggest that SMGs on average contain AGN, although the absence of other AGN tracers in most cases (e.g., Alexander et al. 2005; Pope et al. 2008; Laird et al. 2010; Wang et al. 2013) indicates that these are unlikely to dominate the energetics.

3.6. Metallicity from the [O iii]52 $\mu {\rm{m}}$ and [N iii]57 $\mu {\rm{m}}$ Emission

As shown by Nagao et al. (2011), the [O iii]52/[N iii]57 flux ratio can be a tracer of gas-phase metallicity (see also M. Pereira-Santaella et al. 2017, in preparation). In Figure 5, we compare the measured [O iii]52/[N iii]57 ratio for the stack and for NGP.NA.144, which is the only individual SMG with a detection in at least one of the relevant transitions.

The models in Figure 5 are from Nagao et al. (2011), and show the variation of [O iii]52/[N iii]57 with metallicity for different ionization parameters ( ${\mathrm{log}}_{10}(U)=-2.5$ to −1.5; dimensionless) and gas densities ( ${\mathrm{log}}_{10}(n/{\mathrm{cm}}^{-3})=1$ to 3), compared to the measurements for SMGs. These models are generated with cloudy (Ferland et al. 1998) and include PDRs and H ii regions, although Nagao et al. (2011) show that the [O iii]52 and [N iii]57 emission is mostly from the H ii regions. The range of densities investigated by Nagao et al. (2011) is consistent with the values that we find for the PDRs in average SMGs (Section 3.7). The ionization parameter (U) used by Nagao et al. (2011) to trace the strength of the ionizing source is defined as the ratio of hydrogen-ionizing photons to total hydrogen density. They consider values of U that are typical of H ii regions—the main sources of [O iii]52 and [N iii]57 emission—and are therefore valid for our SMGs. These model U cannot be directly compared with the G₀ from our PDR results (Section 3.7) because U is dependant on total hydrogen density, whereas the hydrogen in PDRs is primarily atomic.²⁶

Metallicity measurements from [O iii]/[N iii] are not expected to be significantly affected by different optical thickness of the [O iii]52 and [N iii]57 lines, because the wavelength difference is small. [O iii]52 and [N iii]57 also have similar filling factors, and the ratio is not affected by differential magnification because the two lines have similar ionization parameters and critical densities to each other. Therefore, they are emitted from approximately the same region of the galaxies. There is unlikely to be a major effect from any weak AGN component, because neither line is significantly boosted by AGN emission. However, the presence of AGN could enhance the ionization parameter, which, as can be seen from Figure 5, would serve to increase the metallicity for a given observed [O iii]52/[N iii]57. However, if AGN emission is dominant (unlikely for SMGs; e.g. Section 3.5), the [O iii]/[N iii] ratio is no longer a good metallicity tracer because the models are for H ii regions and not XDRs.

The results from the stacked data (Figure 5) show that the average [O iii]52/[N iii]57 ratio of SMGs is indicative of them containing enriched gas, with average metallicities $Z\gtrsim {Z}_{\odot }$ . The [O iii]52 and [N iii]57 data are insufficient to constrain the metallicity of NGP.NA.144. Previous measurements of the metallicities of SMGs are hard to come by and have large uncertainties. There have been indications that they typically have approximately metallicities ranging from sub-solar to a few times solar—including from the [N ii] $\lambda 6584$ Å/Hα (Swinbank et al. 2004), extrapolations from the mass–metallicity (or mass–metallicity-SFR) relation, and [N ii]205 μm/[C ii]158 μm for one SMG at $z\sim 4.8$ (Nagao et al. 2012)—although these measurements were plagued with uncertainties and systematic effects.

3.7. ISM Density and FUV Radiation Field

Rather than considering each line observation in isolation, more can be learned by examining the emission line ratios in concert with PDR modeling. However, the low signal-to-noise ratios of the individual observations mean that this is only possible for the stacked data—i.e., we can only examine the properties of the average SMG.

We use the PDR models of Kaufman et al. (1999, 2006), accessed via PDR Toolbox ²⁷ (Pound & Wolfire 2008). The model is characterized using a varying gas density (n, in units of the density of hydrogen nuclei) and the strength of the FUV (energies $h\nu =6$ –13 eV) radiation field (G₀, in units of the Habing Field, $1.6\times {10}^{-3}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ ). Figure 6 highlights where the model n and G₀ produce [O i]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [Si ii]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ consistent with the SMG average values measured from the fiducial mean stacks (Section 3.3 and Table 4). [Fe ii]26 is also available in the Kaufman et al. (2006) PDR model, but our non-detection is too shallow to be useful in constraining n and G₀, so it is not included in Figure 6 or the following discussions.

**Figure 6.** Left: contours showing constraints on the gas density (n) and FUV field strength (G₀ in units of the Habing Field) for the Kaufman et al. (1999, 2006) PDR model from our mean stacked spectra. Thick black lines represent the average of the different line ratios (as labeled), and the colored regions represent the $\pm 1\sigma$ measurement uncertainties. Kaufman et al. (2006) includes two models for [Si ii]34 with metallicity equal to the solar neighborhood value (Z = 1) and for three times solar (Z = 3), with both shown. There is also the possibility that the [O i]63 emission is self-absorbed, so we include vectors on the [O i]/FIR and [O i]/[C ii] contours to show how the positions of the contours would change if the average intrinsic [O i]63 strength is twice that measured (or for [O i]/[C ii] the measured [C ii]158 is twice that from the PDR alone; Section 3.7.1). The dotted blue line labeled "[O i]/[C ii]&FF" shows where the [O i]/[C ii] contour would lie if the measured value were corrected according to the relative filling factors of [O i] and [C ii] from M82 (Section 3.7.1). The green shaded region represents the region of n–G₀ space that we determine is most representative of the average SMG, once the additional considerations described in Section 3.7.1 are taken into account (Section 3.7.2). Right: comparison of the gas density and FUV field strength for average SMGs (derived here; green shaded) compared with local star-forming galaxies (Malhotra et al. 2001) and existing SMG measurements. The existing archival SMG measurements are mostly derived from only [C ii]158 and CO data (Sturm et al. 2010; Cox et al. 2011; Danielson et al. 2011; Valtchanov et al. 2011; Alaghband-Zadeh et al. 2013; Huynh et al. 2014; Rawle et al. 2014) and have typical uncertainties of 0.2–0.5 in both ${\mathrm{log}}_{10}(n/{\mathrm{cm}}^{-3})$ and ${\mathrm{log}}_{10}({G}_{0})$ . The right-hand axis shows approximate sources sizes for a SMG with ${L}_{\mathrm{IR}}={10}^{13}$ ${L}_{\odot }$ at different G₀ values (Section 3.7.2). Contours highlight the regions that typical ULIRGs inhabit, from the GOALS (T. Díaz-Santos et al. 2017, in preparation) and HERUS samples (Farrah et al. 2013), as well as GOALS star-forming galaxies (T. Díaz-Santos et al. 2017, in preparation). We also show the areas that are found to be dominated by local starbursts, local spiral galaxies and Galactic molecular clouds, and Galactic OB star-forming regions on the basis of [C ii]158/CO(J = 1 → 0) measurements (Stacey et al. 1991). The density and FUV field strength from our fine-structure line data are similar to individual local star-forming galaxies, but lower in density than the majority of individually measured SMGs (discussed in Section 3.7.2). The comparison with GOALS and HERUS shows that we do not have sufficient data to robustly distinguish whether the SMGs have internal conditions more similar to local ULIRGs or star-forming galaxies.
Download figure:
Standard image High-resolution image

$\pm 1\sigma $ — **Figure 6.** Left: contours showing constraints on the gas density (n) and FUV field strength (G₀ in units of the Habing Field) for the Kaufman et al. (1999, 2006) PDR model from our mean stacked spectra. Thick black lines represent the average of the different line ratios (as labeled), and the colored regions represent the $\pm 1\sigma$ measurement uncertainties. Kaufman et al. (2006) includes two models for [Si ii]34 with metallicity equal to the solar neighborhood value (Z = 1) and for three times solar (Z = 3), with both shown. There is also the possibility that the [O i]63 emission is self-absorbed, so we include vectors on the [O i]/FIR and [O i]/[C ii] contours to show how the positions of the contours would change if the average intrinsic [O i]63 strength is twice that measured (or for [O i]/[C ii] the measured [C ii]158 is twice that from the PDR alone; Section 3.7.1). The dotted blue line labeled "[O i]/[C ii]&FF" shows where the [O i]/[C ii] contour would lie if the measured value were corrected according to the relative filling factors of [O i] and [C ii] from M82 (Section 3.7.1). The green shaded region represents the region of n–G₀ space that we determine is most representative of the average SMG, once the additional considerations described in Section 3.7.1 are taken into account (Section 3.7.2). Right: comparison of the gas density and FUV field strength for average SMGs (derived here; green shaded) compared with local star-forming galaxies (Malhotra et al. 2001) and existing SMG measurements. The existing archival SMG measurements are mostly derived from only [C ii]158 and CO data (Sturm et al. 2010; Cox et al. 2011; Danielson et al. 2011; Valtchanov et al. 2011; Alaghband-Zadeh et al. 2013; Huynh et al. 2014; Rawle et al. 2014) and have typical uncertainties of 0.2–0.5 in both ${\mathrm{log}}_{10}(n/{\mathrm{cm}}^{-3})$ and ${\mathrm{log}}_{10}({G}_{0})$ . The right-hand axis shows approximate sources sizes for a SMG with ${L}_{\mathrm{IR}}={10}^{13}$ ${L}_{\odot }$ at different G₀ values (Section 3.7.2). Contours highlight the regions that typical ULIRGs inhabit, from the GOALS (T. Díaz-Santos et al. 2017, in preparation) and HERUS samples (Farrah et al. 2013), as well as GOALS star-forming galaxies (T. Díaz-Santos et al. 2017, in preparation). We also show the areas that are found to be dominated by local starbursts, local spiral galaxies and Galactic molecular clouds, and Galactic OB star-forming regions on the basis of [C ii]158/CO(J = 1 → 0) measurements (Stacey et al. 1991). The density and FUV field strength from our fine-structure line data are similar to individual local star-forming galaxies, but lower in density than the majority of individually measured SMGs (discussed in Section 3.7.2). The comparison with GOALS and HERUS shows that we do not have sufficient data to robustly distinguish whether the SMGs have internal conditions more similar to local ULIRGs or star-forming galaxies.
Download figure:
Standard image High-resolution image

The [O i]63 and [Si ii]34 data alone cannot constrain the conditions of the PDRs in SMGs, so we also include archival [C ii]158 measurements from Gullberg et al. (2015) and George (2015) for gravitationally lensed SMGs. Many of the sources in these papers were included in our PACS stacks. The mean [C ii]/FIR ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [O i]/[C ii] are shown in Figure 6, where both the [C ii]158 and [O i]63 fluxes are scaled by the FIR (42.5–122.5 μm) of the sources measured, so as to remove any luminosity effects.

To interpret ISM conditions via spectroscopy and PDR models, it is typical to identify the regions of n and G₀ space where the constraints from different line measurements overlap (e.g., Kaufman et al. 1999, 2006; Pound & Wolfire 2008). It is important to note that the PDR models assume that all the measured fluxes are being emitted from the same spatial region. However, for our SMGs, PACS cannot resolve different regions, and thus each line measurement is an aggregate over the whole galaxy. Because different PDRs within a galaxy may have different properties, the different line observations may, therefore, be differently weighted toward different regions. Furthermore, we are investigating stacked data—i.e., average line strengths over many SMGs—and therefore emission from several galaxies, which may also have intrinsic spread in their properties (Section 3.4). Therefore, our conclusions are averages, with some natural weighting toward more line-luminous regions and galaxies. In Section 3.7.1, we discuss further considerations; our final interpretation of the PDR parameters are discussed in Section 3.7.2.

3.7.1. Additional Considerations

There are several factors that affect the interpretation of Figure 6, which we now discuss. First, a small fraction of local galaxies exhibit self-absorption in the [O i]63 line (e.g., Fischer et al. 1997, 1999; Genzel & Cesarsky 2000; Farrah et al. 2013; Rosenberg et al. 2015), and it is possible that the [O i]63 emission of SMGs may be self-absorbed because our data are too shallow to directly measure this via the shape of the line. Therefore, in Figure 6, we demonstrate that a factor of two increase in the [O i]63 flux from the measured value would have only a minor effect on the positioning of the [O i]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [O i]/[C ii] contours; both are within $2\sigma$ of the directly measured values. This factor of two demonstrates the approximate maximum change in Figure 6 likely from [O i]63 self-absorption.

Second, as discussed in Section 3.3, there is a difference in the [O i]63 emission as measured from the (fiducial) mean stack, and the ${FIR}$ -weighted mean stack. If we use the limit on [O i]63 from the weighted mean stack (instead of the fiducial mean stack; Figure 6), the acceptable [O i]/ $\mathrm{FIR}$ ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [O i]/[C ii] regions would both be substantially larger than the fiducial result, as shown in Figure 7. Overall, considering the weighted stacked fluxes expands the acceptable G₀ range (before accounting for size arguments; Section 3.7.2), but has minimal effect on the inferred PDR density.

**Figure 7.** Contours showing the effect of the constraints on the gas density (n) and FUV field strength (G₀ in units of the Habing Field) for the Kaufman et al. (1999, 2006) PDR model, if the $3\sigma$ limit on [O i]63 emission from the FIR-weighted stacked is used rather than the fiducial mean stack (shown in Figure 6). Colors are the same as Figure 6, with shading showing acceptable regions based on each stacked measurement, with the addition of black arrows to highlight that both the [O i]/ $\mathrm{FIR}$ and [O i]/[C ii] are $3\sigma$ upper limits.
Download figure:
Standard image High-resolution image

$3\sigma $ — **Figure 7.** Contours showing the effect of the constraints on the gas density (n) and FUV field strength (G₀ in units of the Habing Field) for the Kaufman et al. (1999, 2006) PDR model, if the $3\sigma$ limit on [O i]63 emission from the FIR-weighted stacked is used rather than the fiducial mean stack (shown in Figure 6). Colors are the same as Figure 6, with shading showing acceptable regions based on each stacked measurement, with the addition of black arrows to highlight that both the [O i]/ $\mathrm{FIR}$ and [O i]/[C ii] are $3\sigma$ upper limits.
Download figure:
Standard image High-resolution image

In addition, the Kaufman et al. (2006) model includes two sets of [Si ii]34 data, for two different metallicities—labelled Z = 1 (where the gas phase metallicities are those in the solar neighborhood) and Z = 3 (where all elements are three times more abundant) in Figure 6. As discussed in Section 3.6, SMGs are likely to have $Z\gtrsim {Z}_{\odot }$ , although we cannot distinguish between $Z={Z}_{\odot }$ and $Z=3{Z}_{\odot }$ with current data. Due to the uncertainties, we include the constraints for both model metallicities in Figure 6.

For the main contours on Figure 6, the relative filling factors of the various line species in the PACS beams are considered to be equal, i.e., we have not applied any corrections for filling factors. Such corrections are expected to have the biggest effect on the [O i]/[C ii] ratio. To make corrections due to the relative [O i]63 and [C ii]158 filling factors, we would need to know the relative sizes of the regions that dominate those emission lines. The large PACS (for [O i]63) and APEX or SPIRE (for [C ii]158) beam sizes preclude directly measuring the extent of the emission. Instead, to gain some insight into the possible size of this effect, we consider the local starburst M82, where the extents of the [O i]63 and [C ii]158 emission regions can be directly measured, leading to a required correction factor of 0.112 on the [C ii]158 flux, i.e., the [O i]/[C ii] is increased by a factor of $1/0.112=8.9$ (e.g., Stacey et al. 1991; Lord et al. 1996; Kaufman et al. 1999; Contursi et al. 2013). Under these circumstances, the low density ( $n\lesssim {10}^{3}\,{\mathrm{cm}}^{-3}$ ) end of the [O i]/[C ii] contour on Figure 6 would be shifted to higher G₀ (demonstrated with the dotted line on Figure 6), with the $1\sigma$ uncertainty region encompassing ${G}_{0}={10}^{1-6.5}$ , ${10}^{2-5.5}$ , and ${10}^{1.8-4.5}$ for $n={10}^{\mathrm{1,2,3}}\,{\mathrm{cm}}^{-3}$ , respectively. The M82 filling factor correction is substantial, and thus the correction to the [O i]/[C ii] ratio for the average SMG is likely to be lower than this. Therefore, the correction explored here demonstrates an approximate upper boundary to the size of the effect.

Another related consideration is that some of the line emission may originate from H ii regions (or other gas), rather than PDRs. This is most likely to affect the [C ii]158 flux due to the critical densities and ionization parameters of the different transitions studied here. Therefore, the non-PDR [C ii]158 emission should be subtracted from the observed [C ii]158 intensity prior to using it to analyze the PDR conditions. This is typically done by using multi-phase modeling (e.g., cloudy; Ferland et al. 1998), or using other transitions (such as [N ii]122 μm) to determine the contribution from H ii regions. However, there are few observations of [N ii]122 in SMGs (Ferkinhoff et al. 2011; Combes et al. 2012; Decarli et al. 2012; Nagao et al. 2012), and the complexity of cloudy modeling coupled with the limitations in our data means that cloudy analysis will not improve the uncertainties in our analysis. We instead investigate the effect of some of the [C ii]158 emission coming from non-PDR gas qualitatively. Note that, if some of the observed [C ii]158 flux is not from the PDRs, then the correct values of [C ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [O i]/[C ii] to use in Figure 6 would be decreased and increased, respectively. The minor effect of a factor of two increase in [O i]/[C ii] is shown on Figure 6. Note that, due to the critical density and ionization potential of [C ii]158, it is unlikely to be dominated by non-PDR emission, i.e., the factor of two considered is approximately the upper limit of any non-PDR correction required. A decrease of [C ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ also has a small effect on Figure 6, shifting the contours to slightly higher G₀ and n, but remaining within the current $1\sigma$ uncertainty area—although the updated error region marginally overlaps with the [O i]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ uncertainties. Thus, different non-PDR emission is unlikely to have a substantial effect on the interpretation of the [C ii]158 and [O i]63 fluxes in Figure 6.

In some cases, [Si ii]34 can be boosted by AGN emission (e.g., Dale et al. 2006, 2009). If some of the observed [Si ii]34 flux is from AGN, then removing this contribution would move the [Si ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ contours on Figure 6 outwards, i.e., similarly to increasing metallicity. As discussed in Section 3.5, there is evidence from the [Si ii]/[S iii] ratio that there is some AGN contribution in the stacked [Si ii]34 data, which may explain why the contours for n and G₀ for our observed [Si ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ (particularly for $Z={Z}_{\odot }$ ) are offset from those derived from [O i]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ and [O i]/[C ii].

It is possible that the different transitions studied may have different optical depths. Although all the lines in Figure 6 are in the IR, they cover a significant range in rest-frame wavelength (34 μm for [Si ii], to 158 μm for [C ii]), and we are examining some of the dustiest galaxies in the universe. Similarly to the filling factors, this effect is most likely to affect the positioning of the [O i]/[C ii] contours, and will shift them in a similar manner (i.e., toward higher G₀ for a given n), due to the shorter wavelength [O i]63 being more strongly affected.

Finally, $\sim 60 \%$ of the galaxies included in the stacks are known to be gravitationally lensed, and we have so far assumed that the lensing amplification is equal in all components of emission. In fact, because these are composed of galaxy-galaxy lenses, differential magnification, caused by different regions of the background galaxy being amplified by different amounts, is possible. If differential magnification is a random effect, then it is more likely to affect analyses of individual galaxies (e.g., in Figure 3) than the average values examined in Figure 6 and for the PDR modeling, where the effects will be minimized due to averaging many sources. However, there may be systematic effects in regions emitting the majority of the different line species, resulting in them being differentially amplified. This is likely to be an important effect, due to biases in the identification of lensed SMGs—a crucial step in the selection of many of the PACS targets.

The spatial resolution of the spectroscopy is insufficient to resolve and model the lensing for each transition individually, so even if we had attained detections for several galaxies individually, we would be unable to determine the differential magnification on a case-by-case basis. Instead, we consider the simulations of Serjeant (2012), who investigated systematic effects in the differential magnification of a simple dusty galaxy model with a variety of foreground galaxy lenses and alignments. Serjeant (2012) found no systematic differential magnification effects in [C ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ for lensed SMGs, and claimed that that systematic effects are similarly unlikely to be found in [O i]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ because [O i]63 and [C ii]158 are observed to be co-spatial in M82 on small scales (Sturm et al. 2010). In that scenario, we would also not expect differential amplification effects in [O i]/[C ii], or [Si ii]/ $\mathrm{FIR}{}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ , which traces the same environments. However, we consider how the contours on Figure 6 would change if half the [C ii]158 emission is instead from a more extended region (e.g., an H ii region) than the PDRs. In this case, the [C ii]158 flux from the smaller PDR will typically be more highly magnified than the extended [C ii]158, and thus the overall effect will be to minimize the fraction of detected [C ii]158 from non-PDR regions, somewhat canceling out the effect of non-PDR [C ii]158 emission on the line ratios. In this case, the [Si ii]34 and [O i]63 emitting gas will be situated in the PDRs along with the [C ii]158, and thus be co-spatial on all but the smallest scales, making them unlikely to be strongly systematically affected by differential magnification. We reiterate that this discussion of the systematic effects of differential amplification is applicable only to the average (i.e., stacked) values of the sample, and that individual galaxies may have quite substantial differential magnification effects.

3.7.2. Inferred PDR Parameters

Based on the data shown in Figure 6 and the discussion in Section 3.7.1, we determine that the average SMG has $n\sim {10}^{1-3.5}\,{\mathrm{cm}}^{-3}$ and ${G}_{0}\sim {10}^{2.2-4.5}$ . In Table 4, we use these values and PDR Toolbox to predict the average strength of the [C i]370 μm, [C i]609 μm, [O i]145 μm, and [Fe ii]26 μm lines from the PDRs in SMGs.

The $n\sim {10}^{4}\,{\mathrm{cm}}^{-3}$ and ${G}_{0}\sim {10}^{-0.2}$ solution where the [O i]/FIR ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ , [C ii]/FIR ${}_{(42.5\mbox{--}122.5\mu {\rm{m}})}$ , and [O i]/[C ii] contours are coincident (Figure 6) is excluded because such a low G₀ in SMGs with (intrinsic) ${L}_{\mathrm{IR}}\sim {10}^{12.5-13}$ ${L}_{\odot }$ would require source sizes of hundreds of kpc, which is clearly unphysical.²⁸ For ${G}_{0}\sim {10}^{2.2-4.5}$ and ${L}_{\mathrm{IR}}\sim {10}^{13}$ ${L}_{\odot }$ , the source sizes are expected to be ∼1–10 kpc; lower G₀ values are excluded, as they would require excessively large sources. Measurements show that SMGs have typical total sizes of ∼0.5–10 kpc (Tacconi et al. 2006; Younger et al. 2008; Swinbank et al. 2010; Ivison et al. 2011; Riechers et al. 2011b; Bussmann et al. 2013; Calanog et al. 2014; Ikarashi et al. 2015), usually smaller in the far-IR continuum (dust) than the rest-frame optical (stars) or radio emission (e.g., Simpson et al. 2015). Higher G₀ values are acceptable if the weighted stack for the [O i]63 line is used, although for ${G}_{0}\gtrsim {10}^{5}$ the inferred source sizes would be smaller than typically observed for SMGs, and are thus unlikely.

In the right-hand panel of Figure 6, the values of n and G₀ that we infer for average SMGs are compared with the regions of n–G₀ space typically populated by local ULIRGs and star-forming galaxies, as determined by similar PDR modeling from HERUS (Farrah et al. 2013) and GOALS (T. Díaz-Santos et al. 2017, in preparation). The HERUS and GOALS ULIRGs include galaxies with ${L}_{\mathrm{IR}}\geqslant {10}^{12}$ ${L}_{\odot }$ (with $\gt 95 \%$ star formation rather than AGN-dominated); the GOALS star-forming sample consists of the galaxies with $\gt 50 \%$ contribution to the bolometric luminosity from star formation, of which $\sim 90 \%$ are LIRGs and the remainder are ULIRGs. Also highlighted in Figure 6 are the regions found to be preferentially occupied by local starbursts, spiral galaxies, molecular clouds, and galactic OB star-forming regions, as determined by [C ii]158/CO(J = 1 → 0) ratios (Stacey et al. 1991). In addition, we also show measurements of local star-forming galaxies (Malhotra et al. 2001) and individual SMGs with existing measurements of n and G₀, typically from CO and [C ii]158 lines (Sturm et al. 2010; Cox et al. 2011; Danielson et al. 2011; Valtchanov et al. 2011; Alaghband-Zadeh et al. 2013; Huynh et al. 2014; Rawle et al. 2014).

It can be seen in Figure 6 that full PDR modeling (from the HERUS and GOALS results; Farrah et al. 2013; T. Díaz-Santos et al. 2017, in preparation) results in correlated G₀ and n. The local samples also show that local ULIRGs typically have higher G₀ than the mostly LIRG and sub-LIRG star-forming galaxies (although there is substantial overlap between the two populations). This is in contrast with the [C ii]158/CO(J = 1 → 0) measurements from Stacey et al. (1991), which suggested that G₀ can efficiently distinguish between starburst and spiral galaxy star formation. The observations of individual local, mostly (sub-)LIRG, star-forming galaxies (primarily based on [C ii]158 and [O i]63 data; Malhotra et al. 2001) extend to higher n and G₀ than the GOALS results, although they broadly follow the same trend. The PDR density and FUV field strength from our stacked line measurements of SMGs align with the regions for both local ULIRGs and local star-forming galaxies (Figure 6). However, the current data are unable to distinguish between probable merger triggering and secular evolution of SMGs.

Figure 6 also shows that the n values derived from our fine-structure line spectroscopy are lower than most archival measurements of individual SMGs (both lensed and unlensed). It is possible that the apparent difference between our stacking results and individual archival SMG studies are due to bias in the selection of individually analyzed SMGs or uncertainties in our analysis, such as if [O i]63 self-absorption is significantly greater than the factor of two that we have investigated. However, it could also be due to systematic uncertainties in the archival analyses, because most of the archival measurements are inferred from only two spectral features, CO and [C ii]158, and thus n is primarily constrained by CO. These existing SMG studies typically use high-J CO observations, which are converted to CO(J = 1 → 0) fluxes using standard ratios. However, in SMGs, CO(J = 1 → 0) and high-J emission regions are often not co-located, with higher-J lines tracing warmer, more compact dust than CO(J = 1 → 0)(e.g., Ivison et al. 2011; Riechers et al. 2011a; Spilker et al. 2015). This means that, for SMGs, CO(J = 1 → 0) fluxes estimated from higher-J measurements may not trace the same region or spatial scales as the [C ii]158 and PDR models, which could bias the derived n. Bisbas et al. (2014), for example, demonstrate a similar effect, showing that different CO transitions probe to different cloud depths in the PDRs in NGC 4030.

4. Conclusions

We have presented Herschel-PACS spectroscopy and photometry of a sample of 13 gravitationally lensed SMGs at z = 1.03–3.27, targeting the [O iv]26, [Fe ii]26, [S iii]33, [Si ii]34, [O iii]52, [N iii]57, and [O i]63 fine-structure lines, and the ${{\rm{H}}}_{2}$ S(0) and ${{\rm{H}}}_{2}$ S(1) hydrogen rotational lines. We detected only two lines at $\geqslant 3\sigma$ significance ([O iii]52 in NGP.NA.144, and [O iv]26 in HXMM01).

To supplement our data, we identified 32 additional SMGs that also have Herschel-PACS spectroscopy of the targeted lines, and we stacked these archival spectra with those from our 13 originally targeted sources. The stacked spectra include eight (for the two hydrogen lines) to 37 (for [O iii]52) SMGs, resulting in average spectra of SMGs with improvements of up to a factor of ∼6 in the nominal noise level. We detected [O i]63, [Si ii]34, and [N iii]57 in the stacks, with line strengths relative to the far-IR continuum of $(0.36\pm 0.12)\times {10}^{-3}$ , $(0.84\pm 0.17)\times {10}^{-3}$ , and $(0.27\pm 0.10)\times {10}^{-3}$ , respectively.

Based on the [O iii]52/[N iii]57 line ratios we determined that SMGs are typically enriched galaxies, with gas-phase metallicities $\gtrsim {Z}_{\odot }$ . The stacked [Si ii]34/[S iii]33 ratio indicates that there is some LINER and/or AGN contribution to the IR spectra, although the absence of strong [O iv]26 emission and other AGN tracers suggests that these are unlikely to dominate the energetics of typical SMGs.

The ratio of the [O i]63 flux to the far-IR continuum flux in the stacked data is $(0.36\pm 0.12)\times {10}^{-3}$ , significantly lower than the roughly 2 × 10⁻³ observed in local sub-LIRGs, but consistent with local ULIRGs. We used PDR Toolbox to model the stacked [O i]63 and [Si ii]34 data and also included average [C ii]158 measurements from Gullberg et al. (2015) and George (2015). The model indicates that, on average, the PDRs in SMGs have gas densities, n, of ${10}^{1-3}\,{\mathrm{cm}}^{-3}$ and FUV field strengths, ${G}_{0}={10}^{2.2-4.5}$ . These values are consistent with both measurements of local ULIRGs, and most LIRG and sub-LIRG star-forming galaxies. Additional IR data are required to further constrain the PDR models and determine whether the star formation in high-redshift SMGs is more similar to local sub-LIRGs than local ULIRGs. The derived n is lower than most measurements of individual SMGs from [C ii]158 and CO data, which may be due to the previous widespread use of high-J CO transitions and the uncertainties converting these to CO(J = 1 → 0) luminosities.

The authors thank Lee Armus, Jeronimo Bernard-Salas, Kristen Coppin, Tanio Díaz-Santos, Jonatan Selsing, Ian Smail, Mark Swinbank, and Zhi-Yu Zhang for helpful discussions, and an anonymous referee for constructive and clarifying comments. We are grateful to Javier Graciá-Carpio, Eckhard Sturm, and the SHINING collaboration for providing a compilation of spectral observations and IR emission measurements used for comparison with our sample, and to Lucia Marchetti, Steve Crawford, Stephen Serjeant, and Andrew Baker for discussion regarding a followup South African Large Telescope program targeting lens redshifts (L. Marchetti et al. 2017, in preparation).

J.L.W. is supported by a European Union COFUND/Durham Junior Research Fellowship under EU grant agreement number 267209, and acknowledges additional support from STFC (ST/L00075X/1). Support for this work was provided by NASA through an award issued by JPL/Caltech. A.C. and J.L.W. acknowledge support from NSF AST-1313319 and NASA NNX16AF38G. H.D. acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2014 Ramón y Cajal program MINECO RYC-2014-15686. L.D., S.J.M. and R.J.I. acknowledge support from European Research Council Advanced Investigator Grant COSMICISM, 321302; S.J.M. and L.D. are also supported by the European Research Council Consolidator Grant CosmicDust (ERC-2014-CoG-647939, PI: H.L. Gomez). A.V. acknowledges support from the Leverhulme Trust through a Research Fellowship. The Dark Cosmology Centre is funded by the Danish National Research Foundation. J.L.W., A.C., H.D., and D.R. thank the Aspen Center for Physics for hospitality. This work is supported in part by the NSF under Grant Numbers PHY-1066293 and AST-1055919. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe."

This research has made use of data from the HerMES project (http://hermes.sussex.ac.uk). HerMES is a Herschel Key Programme utilizing Guaranteed Time from the SPIRE instrument team, ESAC scientists, and a mission scientist. The data presented in this paper will be released through the HerMES Database in Marseille, HeDaM (http://hedam.oamp.fr/HerMES). The Herschel-ATLAS (http://www.h-atlas.org) is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia, with important participation from NASA. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, and Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). HCSS/HSpot/HIPE are joint developments of the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia. This research has made use of NASA's Astrophysics Data System.

Facilities: Herschel (PACS), Herschel (SPIRE).

Appendix A: Archival Data Included in These Analyses

The archival sources and PACS observations that are included in the stacking analyses are listed in Table 5.

Appendix B: 1D Spectra

In Figure 8 (Figures 8.1–8.13), we show the individual spectral observations for each lensed SMG observed in OT2_jwardlow_1.

The Interstellar Medium in High-redshift Submillimeter Galaxies as Probed by Infrared Spectroscopy^{^∗}

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations and Data Reduction

2.1. Targeted Sample Selection

2.2. Herschel-PACS Spectroscopy

2.3. Herschel-PACS Photometry

2.4. Archival Sample and Data

3. Analysis and Discussion

3.1. Far-Infrared SED Fits

3.2. Individual Emission Line Measurements

3.3. Stacked Spectra

3.4. Individual Line Strengths

3.4.1. [O i]63 $\mu {\rm{m}}$

3.4.2. [S iii]33 $\mu {\rm{m}}$

3.4.3. [Si ii]34 $\mu {\rm{m}}$

3.4.4. [O iii]52 $\mu {\rm{m}}$

3.4.5. [N iii]57 $\mu {\rm{m}}$

3.4.6. [O iv]26 $\mu {\rm{m}}$

3.4.7. [Fe ii]26 $\mu {\rm{m}}$

3.4.8. H₂ S(0) and H₂ S(1)

3.5. AGN Contribution Traced by [Si ii]34 $\mu {\rm{m}}$ and [S iii]33 $\mu {\rm{m}}$ Emission

3.6. Metallicity from the [O iii]52 $\mu {\rm{m}}$ and [N iii]57 $\mu {\rm{m}}$ Emission

3.7. ISM Density and FUV Radiation Field

3.7.1. Additional Considerations

3.7.2. Inferred PDR Parameters

4. Conclusions

Appendix A: Archival Data Included in These Analyses

Appendix B: 1D Spectra

Footnotes

The Interstellar Medium in High-redshift Submillimeter Galaxies as Probed by Infrared Spectroscopy∗

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations and Data Reduction

2.1. Targeted Sample Selection

2.2. Herschel-PACS Spectroscopy

2.3. Herschel-PACS Photometry

2.4. Archival Sample and Data

3. Analysis and Discussion

3.1. Far-Infrared SED Fits

3.2. Individual Emission Line Measurements

3.3. Stacked Spectra

3.4. Individual Line Strengths

3.4.1. [O i]63 \mu {\rm{m}}

3.4.2. [S iii]33 \mu {\rm{m}}

3.4.3. [Si ii]34 \mu {\rm{m}}

3.4.4. [O iii]52 \mu {\rm{m}}

3.4.5. [N iii]57 \mu {\rm{m}}

3.4.6. [O iv]26 \mu {\rm{m}}

3.4.7. [Fe ii]26 \mu {\rm{m}}

3.4.8. H2 S(0) and H2 S(1)

3.5. AGN Contribution Traced by [Si ii]34 \mu {\rm{m}} and [S iii]33 \mu {\rm{m}} Emission

3.6. Metallicity from the [O iii]52 \mu {\rm{m}} and [N iii]57 \mu {\rm{m}} Emission

3.7. ISM Density and FUV Radiation Field

3.7.1. Additional Considerations

3.7.2. Inferred PDR Parameters

4. Conclusions

Appendix A: Archival Data Included in These Analyses

Appendix B: 1D Spectra

Footnotes

The Interstellar Medium in High-redshift Submillimeter Galaxies as Probed by Infrared Spectroscopy^{^∗}

3.4.1. [O i]63 $\mu {\rm{m}}$

3.4.2. [S iii]33 $\mu {\rm{m}}$

3.4.3. [Si ii]34 $\mu {\rm{m}}$

3.4.4. [O iii]52 $\mu {\rm{m}}$

3.4.5. [N iii]57 $\mu {\rm{m}}$

3.4.6. [O iv]26 $\mu {\rm{m}}$

3.4.7. [Fe ii]26 $\mu {\rm{m}}$

3.4.8. H₂ S(0) and H₂ S(1)

3.5. AGN Contribution Traced by [Si ii]34 $\mu {\rm{m}}$ and [S iii]33 $\mu {\rm{m}}$ Emission

3.6. Metallicity from the [O iii]52 $\mu {\rm{m}}$ and [N iii]57 $\mu {\rm{m}}$ Emission