The interstellar medium in high-redshift submillimeter galaxies as probed by infrared spectroscopy

TL;DR

This study uses infrared spectroscopy to analyze the interstellar medium in high-redshift submillimeter galaxies, revealing their metallicities, gas densities, and star formation conditions, and comparing them to local ULIRGs and disk galaxies.

Contribution

It provides the first comprehensive spectral analysis of SMGs' ISM conditions at high redshift using Herschel-PACS data, including stacking spectra for improved detection.

Findings

SMGs have gas-phase metallicities ≥ solar.

Average gas densities are ~10^{1-3} cm^{-3}.

FUV field strengths are ~10^{2.2-4.5} Habing units.

Abstract

Submillimeter galaxies (SMGs) at $z\gtrsim1$ are luminous in the far-infrared and have star-formation rates, SFR, of hundreds to thousands of solar masses per year. However, it is unclear whether they are true analogs of local ULIRGs or whether the mode of their star formation is more similar to that in local disk galaxies. We target these questions by using Herschel-PACS to examine the conditions in the interstellar medium (ISM) in far-infrared luminous SMGs at z~1-4. We present 70-160 micron photometry and spectroscopy of the [OIV]26 micron, [FeII]26 micron, [SIII]33 micron, [SiII]34 micron, [OIII]52 micron, [NIII]57 micron, and [OI]63 micron fine-structure lines and the S(0) and S(1) hydrogen rotational lines in 13 lensed SMGs identified by their brightness in early Herschel data. Most of the 13 targets are not individually spectroscopically detected and we instead focus on stacking these spectra with observations of an additional 32 SMGs from the \herschel\ archive -- representing a complete compilation of PACS spectroscopy of SMGs. We detect [OI]63 micron, [SiII]34 micron, and [NIII]57 micron at >3sigma in the stacked spectra, determining that the average strengths of these lines relative to the far-IR continuum are $(0.36\pm0.12)\times10^{-3}$, $(0.84\pm0.17)\times10^{-3}$, and $(0.27\pm0.10)\times10^{-3}$, respectively. Using the [OIII]52/[NIII]57 emission line ratio we show that SMGs have average gas-phase metallicities $\gtrsim Z_{\rm sun}$. By using PDR modelling and combining the new spectral measurements with integrated far-infrared fluxes and existing [CII]158 micron data we show that SMGs have average gas densities, n, of $\sim10^{1-3}{\rm cm^{-3}}$ and FUV field strengths, $G_0\sim10^{2.2-4.5}$ (in Habing units: $1.6\times10^{-3}{\rm erg~cm^{-2}~s^{-1}}$), consistent with both local ULIRGs and lower luminosity star-forming galaxies.

Tables

Table 1. Table 1 Positions, redshifts and the lensing amplifications of the target galaxies

Target	Short Names	$z_{\mathrm{source}}$	$z_{\mathrm{lens}}$	Magnificationa	Referencesb	OBSIDsc
H-ATLAS J142935.3$-$002836	G15v2.19,	1.027	0.218	$9.7\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 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\pm 0.4$	B13, W13, C14, Rp	134225952[1–4], 13422620[39,40],
						1342257689

Table 2. Table 2 PACS 70- and 160-µm photometry and derived far-IR fluxes and luminosities.

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

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	Continuumb	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb
	(10-18Wm-2)	(µm)	(mJy)	(10-18Wm-2)	(µm)	(mJy)	(10-18Wm-2)	(µm)	(mJy)
G15v2.19		128.07	$1044\pm 157$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
G09v1.40		195.29		$\cdots$	$\cdots$	$\cdots$		107.61
G12v2.257	$\cdots$	$\cdots$	$\cdots$		106.84			111.09
NGP.NA.144		202.30			107.21			111.47
NGP.NA.56		208.62			110.55			114.96
HXMM01		208.94			110.72			115.13
G09v1.124	$\cdots$	$\cdots$	$\cdots$		114.17			118.72
G15v2.235	$\cdots$	$\cdots$	$\cdots$		116.48			121.12
G09v1.326	$\cdots$	$\cdots$	$\cdots$		119.96		$\cdots$	$\cdots$	$\cdots$
NGP.NB.78	$\cdots$	$\cdots$	$\cdots$		137.64	$128\pm 91$		143.12	$151\pm 84$
G12v2.43	$\cdots$	$\cdots$	$\cdots$		138.21	$103\pm 82$		143.71	$119\pm 76$
G12v2.30	$\cdots$	$\cdots$	$\cdots$		142.60			148.27
HBootes01	$\cdots$	$\cdots$	$\cdots$		143.10			148.80
Mean Stackc	$1.0\pm 0.3$	$63.16$	…		$\cdots$	…	$2.8\pm 0.4$	$34.83$	…
	[O III]52µm			[N III]57µm
Name	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb
	(10-18Wm-2)	(µm)	(mJy)	(10-18Wm-2)	(µm)	(mJy)
G15v2.19	$\cdots$	$\cdots$	$\cdots$		116.19	$913\pm 132$
G09v1.40		160.14			177.18	$220\pm 114$
G12v2.257		165.33		$\cdots$	$\cdots$	$\cdots$
NGP.NA.144	$10.6\pm 3.2$	165.86			183.54
NGP.NA.56		171.08	$241\pm 137$		189.27	$342\pm 187$
HXMM01		171.34	$129\pm 42$	$\cdots$	$\cdots$	$\cdots$
G09v1.124		176.67		$\cdots$	$\cdots$	$\cdots$
G15v2.235		180.25		$\cdots$	$\cdots$	$\cdots$
G09v1.326		185.64		$\cdots$	$\cdots$	$\cdots$
NGP.NB.78	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
G12v2.43	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
G12v2.30	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
HBootes01	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
Mean Stackc		$\cdots$	…	$1.9\pm 0.6$	$57.38$	…
	[O IV]26µmd and [Fe II]26µmd
Name	[O iv] flux	[O iv] ${\lambda }_{\mathrm{obs}}$a	Continuumb	[Fe ii] flux	[Fe ii] ${\lambda }_{\mathrm{obs}}$a
	(10-18Wm-2)	(µm)	(mJy)	(10-18Wm-2)	(µm)
G15v2.19	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
G09v1.40		80.03			80.33
G12v2.257		82.61			82.93
NGP.NA.144		82.90	$226\pm 190$		83.21
NGP.NA.56		85.49			85.81
HXMM01	$7.7\pm 2.5$	85.70			85.94
G09v1.124		88.28			88.62
G15v2.235		90.07			90.41
G09v1.326		92.76			93.12
NGP.NB.78	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
G12v2.43		106.87	$53\pm 26$		107.28
G12v2.30		110.27	$56\pm 40$		110.68
HBootes01	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
Mean Stackc		$\cdots$	$\cdots$		$\cdots$
	H2 S(0)			H2 S(1)
Name	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb	Line flux	Line ${\lambda }_{\mathrm{obs}}$a	Continuumb
	(10-18Wm-2)	(µm)	(mJy)	(10-18Wm-2)	(µm)	(mJy)
G15v2.19		57.20		$\cdots$	$\cdots$	$\cdots$
G09v1.40		87.22		$\cdots$	$\cdots$	$\cdots$
G12v2.257	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
NGP.NA.144	$\cdots$	$\cdots$	$\cdots$		54.55
NGP.NA.56		93.18			56.25
HXMM01	$\cdots$	$\cdots$	$\cdots$		56.33
G09v1.124	$\cdots$	$\cdots$	$\cdots$		58.09
G15v2.235	$\cdots$	$\cdots$	$\cdots$		59.26
G09v1.326	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
NGP.NB.78		116.01	$107\pm 43$	$\cdots$	$\cdots$	$\cdots$
G12v2.43		116.49	$86\pm 40$	$\cdots$	$\cdots$	$\cdots$
G12v2.30		120.18			72.55
HBootes01		120.61		$\cdots$	$\cdots$	$\cdots$
Mean Stackc		$\cdots$	$\cdots$		$\cdots$	$\cdots$

Table 4. Table 4 Relative line strengths for average SMGs.

Ratio	Value
Measurementsa
[O I]63/FIR	$\left(3.6\pm 1.2\right)\times {10}^{-4}$
[S III]33/FIR
[Si II]34/FIR	$\left(8.4\pm 1.7\right)\times {10}^{-4}$
[O III]52/FIR
[N III]57/FIR	$\left(2.7\pm 1.0\right)\times {10}^{-4}$
[O IV]26/FIR
[Fe II]26/FIR
$H_2$ S(0)/FIR
$H_2$ S(1)/FIR
[C II]158/FIR	$\left(1.7\pm 1.1\right)\times {10}^{-3}$
[O I]63/[C II]158	$2.2\pm 1.5$
Predictions from the PDR modelb
[C I]609/FIR	(0.01–20)$\times {10}^{-5}$
[C I]370/FIR	(0.01–60)$\times {10}^{-5}$
[O I]145/FIR	(0.01–20)$\times {10}^{-5}$
[Fe II]26/FIR, $Z=Z_{\odot }$	(0.3–50)$\times {10}^{-7}$
[Fe II]26/FIR, $Z=3Z_{\odot }$	(0.9–800)$\times {10}^{-7}$

Table 5. Table 5 Archival sources and observations included in the stacking analyses.

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

Equations

$$S_{\nu }\propto {\nu }^{\beta }{B}_{\nu }(T_{\mathrm{D}}),$$

$$B_{\nu }(T_{\mathrm{D}})\propto \frac{{\nu }^3}{e^{h\nu /k{T}_{\mathrm{D}}}-1},$$

Full text

{}^{*}$${}^{*}affiliationtext: Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.{\dagger}$${\dagger}affiliationtext: julie.wardlow@durham.ac.uk

1. Introduction

Submillimetre galaxies (SMGs), selected from their high flux densities at submillimetre 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 $\gtrsim 10^{12}$ ${\rm L_{\sun}}$, equivalent to local ultraluminous infrared galaxies (ULIRGs), 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), and SMGs with fluxes down to $\sim 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 $\sim 1000$ ${\rm M_{\sun}yr^{-1}}$) and their gas depletion times suggest that their star formation is episodic and that they are observed in a short-lived (timescales $\sim 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; Menéndez-Delmestre et al., 2013; Hayward et al., 2013; Cowley et al., 2015; Narayanan et al., 2015) and with limited data the 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 inaccesible. 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 H2 S(0) (28µm) and H2 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 in order to probe a range of ISM conditions, in terms of ionization potential and critical density, and 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. Since 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 modelling. Finally, our conclusions are presented in Section 4. Throughout this paper we use $\Lambda$CDM cosmology with $\Omega_{\rm M}=0.27$, $\Omega_{\Lambda}=0.73$ and $H_{0}=71\,{\rm km\,s^{-1}\,Mpc^{-1}}$.

2. Observations and data reduction

In this paper we first analyse PACS observations of sources targeted by our Herschel program OT2_jwardlow_1, which are described in Section 2.1. We later combine these with archival spectroscopy for additional SMGs, which are 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}\geq 100$ mJy; Negrello et al., 2010, 2016; Wardlow et al., 2013; Nayyeri et al., 2016, e.g.). Extensive follow-up programs, including CO spectroscopy (e.g., Frayer et al., 2011; Harris et al., 2012, Riechers et al. in prep.), high-resolution (sub)millimeter and radio interferometry (e.g., Bussmann et al., 2013), high-resolution near-IR imaging (e.g., Wardlow et al., 2013; Negrello et al., 2014a; Calanog 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}\geq 100$ mJy and 70 µm fluxes predicted to be $\geq 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 note that many of the redshifts are from broadband instruments used for line searches, which can have up to $\sim 100~{}{\rm km~{}s^{-1}}$ spectral resolution). PACS spectroscopy of six additional Herschel H-ATLAS and HerMES gravitationally lensed galaxies, and other high-redshift galaxies were observed in a separate program and will be presented in Verma et al. (in prep.), though they are included here in our stacking analyses (see Section 2.4).

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 H2 S(0) and H2 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.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 and in those cases [O IV]26 was observed separately.

The data were reduced using the Herschel Interactive Processing Environment (Ott 2010; hipe) v12.1.0 with version 65.0 of the PACS calibration tree.111We have verified that later versions of hipe do not affect the results by comparing a selection of data reduced with our hipe v12.1.0 script, with v14.0.1 pipeline processed versions of the same observations, and find no significant differences in the reduced spectra. 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 “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 photometry222Due to the enhanced spatial scales from gravitational magnification, approximately half of the targets are marginally resolved by PACS. (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 and 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 and 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 (SNR) threshold to mask the sources during filtering. The final maps are each $\sim 3.5\arcmin\times 7.5\arcmin$ in size with coverage $\geq 90\%$ of the maximum in the central $\sim 0.5\arcmin\times 1\arcmin$ area.

PACS photometry is measured in 18 apertures with radii from 2 to 50″ using the “annularSkyAperturePhotometry” task within hipe, with each measurement 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 $<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.333http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb These total flux densities are presented in Table 2, where we also include PACS 100 µm data from George et al. (in prep.; 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$\sigma$ lower than those presented here, since 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 proposals444www.cosmos.esa.int/web/herschel/observing-overview 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, and 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; Verma et al. in prep.), lensed SPT sources (Vieira et al., 2013), and other SMGs.

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 scaling to the wavelength ranges for $L_{\rm 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, decsribed 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.

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

[TABLE]

where $S_{\nu}$ is the flux density, $\nu$ is frequency and $\beta$ is the power law emissivitity index. $B_{\nu}(T_{\rm D})$ is the Planck function, defined as

[TABLE]

for a dust temperature, $T_{\rm D}$, and where $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_{\rm FIR}$) and 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:

• •

$L_{\rm FIR}$ is the luminosity of the rest-frame SED integrated between 40 and 500 µm, and;

• •

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_{\rm FIR}$ and FIR calculated from the modified blackbody SED fits are listed in Table 3.1 and are used in the analysis in the rest of this paper.

However, since 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_{\rm 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_{\rm 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_{\rm FIR}$ from the Dale & Helou (2002) SED fits that are significantly different to the values from the modified blackbody fits. These are G15v2.19, G12v2.43 and HBoötes01, which are 15% higher and 10% and 10% lower for the Dale & Helou (2002) fits, respectively. Only one galaxy has significantly higher FIR from the Dale & Helou (2002) SEDs than the modified blackbody fits, which is G12v2.257 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 upwards 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_{\rm FIR}$ and FIR from the Dale & Helou (2002) and modified blackbody fits are because $L_{\rm 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_{\rm FIR}$ and the slightly lower normalization of some of the Dale & Helou (2002) fits, having lower $L_{\rm FIR}$ for the Dale & Helou (2002) fits compared with the modified blackbody in some cases is not unexpected.

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, and 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 non-linear Levenberg-Marquardt least-squares minimization. We constrain the fits to have non-negative continua and the velocity offsets of the lines are required to be $\leq\pm 800~{}{\rm 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, and 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 unreliable555PACS Data Reduction Guide for Spectroscopy, §7.7: http://herschel.esac.esa.int/twiki/pub/Public/PacsCalibrationWeb/PDRGspec_HIPE14p2.pdf, and 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*-1* 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.2.