THE HERSCHEL COMPREHENSIVE (U)LIRG EMISSION SURVEY (HERCULES): CO LADDERS, FINE STRUCTURE LINES, AND NEUTRAL GAS COOLING^*

The sample was chosen from the IRAS Revised Bright Galaxy Sample (RBGS), which contains all 629 extragalactic sources with IRAS 60 μm flux density S₆₀ > 5.24 Jy in the (IRAS-covered) sky at Galactic latitudes |b| > 5 (Sanders et al. 2003). From the IRAS RBGS we select a sub-sample applying limits both in S₆₀ and L_IR: at luminosities L_IR > 10¹² L_☉ (ULIRGs), all sources with S₆₀ > 11.65 Jy are included, while at luminosities 10¹¹ L_☉ < L_IR < 10¹² L_☉ (LIRGs), sources with S₆₀ > 16.4 Jy are included. From this flux-limited representative parent sample of 32 targets, we removed three LIRGs for which no ground-based CO data are available, with the exception of ESO 173-G015, IRAS 13120−5453, and MCG+12-02-001. The resulting representative flux-limited sample consists of 21 LIRGs and 8 ULIRGs. The sample covers a factor of 32 in L_IR and contains a range of objects including starburst galaxies, AGNs, and composite sources, and covering also a range of IRAS 60/100 μm ratios. The full list of included galaxies and their respective properties can be found in Table 1. The infrared luminosity and the luminosity distance are from Armus et al. (2009).

Table 1. Sample Properties

1	2	3	4	5	6	7	8	9	10	11
NGC 34	11.49	0.78	0.01962	84.1	330	SB	[O i]63, [O i]145, [C ii]	00h11m0667 −12°06'2613	1342199416	KPOT_pvanderw_1
(IRAS 00085−1223)							194–671 μm	00h11m0653 −12°06'2490	1342199253	KPOT_pvanderw_1
MCG+12–02–001	11.50	1.07	0.01570	69.8	200	SB	[O i]63, [O i]145, [C ii]	00h54m0333 +73°04'5983	1342193211	KPOT_pvanderw_1
(IRAS 00506+7248)							194–671 μm	00h54m0356 +73°05'1038	1342213377	KPOT_pvanderw_1
IC 1623	11.71	1.14	0.02007	85.5	250	SB, AGN	[O i]63, [O i]145, [C ii]	01h07m4659 −17°30'2646	1342212532	KPOT_pvanderw_1
(IRAS 01053−1746)							194–671 μm	01h07m4674 −17°30'2605	1342212314	KPOT_pvanderw_1
NGC 1068	11.40	9.07	0.003793	15.9	300	AGN, SB	[O i]63	02h42m4078 −00°00'4716	1342191153	KPGT_esturm_1K
(IRAS 02401−0013)							[O i]145, [C ii]	02h42m4073 −00°00'4224	1342191154	KPGT_esturm_1K
							194–671 μm	02h42m4092 −00°00'4665	1342213445	KPGT_cwilso01_1
NGC 1365	11.00	4.32	0.00546	17.9	250	Sy1, SB	[O i]63	03h33m3631 −36°08'1661	1342191295	KPGT_esturm_1K
(IRAS 03317−3618)							[O i]145, [C ii]	03h33m3626 −36°08'2433	1342191294	KPGT_esturm_1K
							194–671 μm	03h33m3648 −36°08'1932	1342204020	KPOT_pvanderw_1
NGC 1614	11.65	1.50	0.01594	67.8	220	SB	[O i]63, [O i]145, [C ii]	04h33m5979 −08°34'4419	1342190367	KPOT_pvanderw_1
(IRAS 04315−0840)							194–671 μm	04h33m5985 −08°34'4415	1342192831	KPOT_pvanderw_1
IRAS F05189−2524	12.16	0.60	0.04256	187	300	QSO	[O i]63	05h21m0124 −25°21'4316	1342219441	KPGT_esturm_1K
							[O i]145, [C ii]	05h21m0128 −25°21'4215	1342219442	KPGT_esturm_1K
							194–671 μm	05h21m0142 −25°21'4547	1342192832a	KPOT_pvanderw_1
							194–671 μm	05h21m0142 −25°21'4548	1342192833a	KPOT_pvanderw_1
NGC 2146	11.12	6.97	0.00298	17.5	250	SB	[O i]63, [O i]145, [C ii]	06h18m3553 +78°21'2539	1342193210	KPOT_pvanderw_1
(IRAS 06106+7822)							194–671 μm	06h18m3807 +78°21'2506	1342204025	KPOT_pvanderw_1
NGC 2623	11.60	1.15	0.01851	84.1	400	SB, AGN	[O i]63, [O i]145, [C ii]	08h38m2429 +25°45'1672	1342208904	KPOT_pvanderw_1
(IRAS 08354+2555)							194–671 μm	08h38m2414 +25°45'1734	1342219553	KPOT_pvanderw_1
NGC 3256	11.64	4.61	0.00935	38.9	230	SB	[O i]63	10h27m5161 −43°54'1539	1342210383	KPGT_esturm_1K
(IRAS 10257−4338)							[O i]145, [C ii]	10h27m5145 −43°54'2187	1342210384	KPGT_esturm_1K
							194–671 μm	10h27m5149 −43°54'1600	1342201201	KPOT_pvanderw_1
Arp 299 A	11.93	4.84	0.01030	50.7	325	SB, AGN	[O i]63	11h28m3341 +58°33'4604	1342199421	KPGT_esturm_1K
IC 694							[O i]145	11h28m3341 +58°33'4604	1342232602	OT1_shaileyd_1
(IRAS 11257+5850)							[C ii]	11h28m3341 +58°33'4604	1342208906	KPGT_esturm_1K
							194–671 μm	11h28m3341 +58°33'4604	1342199248	KPOT_pvanderw_1
ESO 320–G030	11.17	1.67	0.01078	41.2	350	SB	[O i]63, [O i]145, [C ii]	11h53m1175 −39°07'5175	1342212227	KPOT_pvanderw_1
(IRAS11506−3851)							194–671 μm	11h53m1152 −39°07'5024	1342210861	KPOT_pvanderw_1
NGC 4418	11.19	1.74	0.007268	36.5	163	Sy2	[O i]63	12h26m5451 −00°52'4077	1342187780	KPGT_esturm_1K
(IRAS 12243−0036)							[O i]145, [C ii]	12h26m5457 −00°52'3693	1342210830	KPGT_esturm_1K
							194–671 μm	12h26m5460 −00°52'3654	1342210848	KPGT_esturm_1K
Mrk 231	12.57	1.48	0.04217	192	200	QSO	[O i]63	12h56m1465 +56°52'2413	1342189280	KPGT_esturm_1K
(IRAS 12540+5708)							[O i]145, [C ii]	12h56m1429 +56°52'2340	1342186811	SDP_esturm_3
							194–671 μm	12h56m1429 +56°52'2673	1342210493	KPOT_pvanderw_1
IRAS13120−5453	12.32	1.94	0.03076	144	400	Sy2, SB	[O i]63	13h15m0628 −55°09'2446	1342214628	KPGT_esturm_1K
							[O i]145, [C ii]	13h15m0617 −55°09'2538	1342214629	KPGT_esturm_1K
							194–671 μm	13h15m0611 −55°09'2321	1342212342	KPOT_pvanderw_1
Arp 193	11.73	8.19	0.02330	110	400	SB, L	[O i]63, [O i]145, [C ii]	13h20m3520 +34°08'2458	1342197801	KPOT_pvanderw_1
(IRAS 13183+3423)							194–671μm	13h20m3535 +34°08'2346	1342209853	KPOT_pvanderw_1
NGC 5135	11.30	0.91	0.01369	60.9	150	Sy2, SB	[O i]63, [O i]145, [C ii]	13h25m4396 −29°50'0174	1342190371	KPOT_pvanderw_1
(IRAS 13229−2934)							194–671 μm	13h25m4391 −29°50'0027	1342212344	KPOT_pvanderw_1
ESO 173–G015	11.38	3.61	0.00974	34	200	SB	[O i]63, [O i]145, [C ii]	13h27m2400 −57°29'2363	1342190368	KPOT_pvanderw_1
(IRAS 13242−5713)							194–671 μm	13h27m2395 −57°29'2289	1342202268	KPOT_pvanderw_1
Mrk 273	12.21	1.05	0.03736	173	520	SB, Sy2	[O i]63	13h44m4209 +55°53'0914	1342207801	KPGT_esturm_1K
(IRAS 13428+5608)							[O i]145, [C ii]	13h44m4182 +55°53'0875	1342207802	KPGT_esturm_1K
							194–671 μm	13h44m4210 +55°53'1050	1342209850	KPOT_pvanderw_1
Zw 049.057	11.35	1.05	0.01300	65.4	200	SB	[O i]63, [O i]145, [C ii]	15h13m1318 +07°13'3071	1342190374	KPOT_pvanderw_1
CGCG 049−057							194–671 μm	15h13m1310 +07°13'2919	1342212346	KPOT_pvanderw_1
(IRAS 15107+0724)
Arp 220	12.28	4.87	0.01813	77	504	SB, AGN	[O i]63, [O i]145	15h34m5722 +23°30'1106	1342191304	KPGT_esturm_1K
(IRAS 15327+2340)							[C ii]	15h34m5721 +23°30'1013	1342191306	KPGT_esturm_1K
							194–671 μm	15h34m5711 +23°30'1126	1342190674	KPGT_cwilso01_1
NGC 6240	11.93	1.10	0.02448	116	500	SB, AGN	[O i]63	16h52m5910 +02°24'0358	1342216622	KPGT_esturm_1K
(IRAS 16504+0228)							[O i]145, [C ii]	16h52m5910 +02°24'0279	1342216623	KPGT_esturm_1K
							194–671 μm	16h52m5901 +02°24'0327	1342214831	KPOT_pvanderw_1
IRAS F17207−0014	12.46	1.56	0.04281	198	620	SB, L	[O i]63	17h23m2183 −00°16'5982	1342229692	KPGT_esturm_1K
							[O i]145, [C ii]	17h23m3184 −00°16'5760	1342229693	KPGT_esturm_1K
							194–671 μm	17h23m2193 −00°17'0110	1342192829	KPOT_pvanderw_1
IRAS F18293−3413	11.88	1.82	0.01817	86	270		[O i]63, [O i]145, [C ii]	18h32m4134 −34°11'3690	1342192112	KPOT_pvanderw_
							194–671 μm	18h32m4117 −34°11'2723	1342192830	KPOT_pvanderw_1
IC 4687/6	11.62	0.84	0.01735	81.9	230	SB	[O i]63	18h13m3994 −57°43'4966	1342239740	OT1_larmus_1
(IRAS 18093−5744)							[C ii]	18h13m3980 −57°43'3571	1342239739	OT1_larmus_1
							194–671 μm	18h13m3950 −57°43'3105	1342192993	KPOT_pvanderw_1
NGC 7469	11.65	1.32	0.01632	70.8	300	Sy1, SB	[O i]63	23h03m1547 +08°52'3705	1342187847	KPGT_esturm_1K
(IRAS 23007+0836)							[O i]145, [C ii]	23h03m1583 +08°52'2852	1342211171	KPGT_esturm_1K
							194–671 μm	23h03m1587 +08°52'2815	1342199252	KPOT_pvanderw_1
NGC 7552	11.11	3.64	0.00537	23.5	180	SB	[O i]63	23h16m1010 −42°34'5389	1342210400	KPGT_esturm_1K
(IRAS 23134−4251)							[O i]145, [C ii]	23h16m1059 −42°35'0573	1342210399	KPGT_esturm_1K
							194–671 μm	23h16m1073 −42°35'0602	1342198428	KPOT_pvanderw_1
NGC 7771	11.40	1.08	0.01427	61.2	250	SB	[O i]63, [O i]145, [C ii]	23h51m2483 +20°06'4233	1342197839	KPOT_pvanderw_1
(IRAS 23488+1949)							194–671 μm	23h51m2472 +20°06'4211	1342212317	KPOT_pvanderw_1
Mrk 331	11.50	0.87	0.01790	79.3	215	SB	[O i]63, [O i]145, [C ii]	23h51m2676 +20°35'0983	1342197840	KPOT_pvanderw_1
(IRAS 23488+2018)							194–671 μm	23h51m2665 +20°35'1042	1342212316	KPOT_pvanderw_1

Notes. Column 1: object name. Column 2: log(L_IR/L_☉) from Armus et al. (2009). Observations use the cosmological parameters H₀ = 70 km s⁻¹ Mpc, Ω_vaccum = 0.72, and Ω_matter = 0.28. Column 3: far infrared flux (FIR) calculated using the IRAS definition (Helou et al. 1985) in 10⁻¹² W m⁻². Column 4: redshift z from NED. Column 5: luminosity distance D_L in Mpc from Armus et al. (2009). Column 6: CO 1–0 full width to half power line width in km s⁻¹. Column 7: galaxy classification from NED SB = starburst, L = LINER, AGN = active galaxy nucleus, Sy1 = Seyfert 1, Sy2 = Seyfert 2, QSO = quasi-stellar object. Column 8: line names. Column 9: pointing coordinates. Column 10: observation ID (OBSID). Column 11: program ID. ^aThe two SPIRE observations of IRAS F05189−2524 were combined using HIPE average (avg) task.

Download table as:  ASCIITypeset images: 1 2

In order to obtain a comprehensive view of the CO emission and the cooling budget of these galaxies, we proposed Herschel/SPIRE spectroscopy (for the CO ladder) and Herschel/PACS spectroscopy (for the [C ii] and [O i] fine structure lines) of the entire sample, unless PACS observations were already observed as part of another program. In addition to the galaxies observed for HerCULES, we have included NGC 4418, NGC 1068, and Arp 220 for completeness. This project was approved as a Key Project on the Herschel Space Observatory, under the name HerCULES (PI: van der Werf). Key elements of HerCULES are:

• 1.  
  a representative flux-limited sample of local LIRGs and ULIRGs;
• 2.  
  comprehensive coverage of the SPIRE spectral range at the highest spectral resolution mode (covering the CO ladder, [C i] and [N ii] fine structure lines, and any other bright features such as H₂O lines);
• 3.  
  comprehensive coverage of the key fine-structure cooling lines [C ii] and [O i] with PACS observations.

Details about the galaxy type and observation ID can be seen in Table 1. We have included observations from other programs (KPGT_esturm_1K, KPGT_cwilso01_1, OT1_larmus_1, OT1_shaileyd_1) to help realize the complete flux-limited sample. The references for these observations are also in Table 1.

Spectra were obtained with the Spectral and Photometric Imaging Receiver and Fourier-Transform Spectrometer (SPIRE-FTS; Griffin et al. 2010) on board the Herschel Space Observatory (Pilbratt et al. 2010) for the full HerCULES sample. The observations were carried out in staring mode with the galaxy nucleus on the central pixel of the detector array, with a beam size varying from 17''–42'' for the CO transitions. The high spectral resolution mode was used with a resolution of 1.2 GHz over the two observing bands. The low frequency focal plane array (Long Wavelength Spectrometer Array, SLW) covers ν = 447–989 GHz (λ = 671–303 μm) and the high frequency focal plane array (Short Wavelength Spectrometer Array, SSW) covers ν = 958–1545 GHz (λ = 313–194 μm), and together they include the CO J = 4–3 to CO J = 13–12 lines. For galaxies with z > 0.03, the rest frequency of the J = 4–3 transition falls short of the SPIRE coverage. All galaxies were observed in the sparse observing mode besides NGC 4418, which was observed in the intermediate mode.

The data were reduced using version 13.0 of the Herschel Interactive Processing Environment (HIPE). Initial processing steps included timeline deglitching, linearity correction, clipping of saturated points, time-domain phase correction, and interferogram baseline subtraction. After a second deglitching step and interferogram phase correction, the interferograms were Fourier transformed, and the thermal emission from instrument and telescope was removed from the resulting spectra. The averaged spectra were flux calibrated as point sources using the calibration tree associated with HIPE 13.0. Following these steps a "dark" spectrum was subtracted, to remove any residual emission from the telescope and the instrument. Since the emission of most of our sources is contained entirely in the central pixel of the detector arrays, a "dark" spectrum was constructed by spectrally smoothing a combination of several off-axis pixels. For extended targets, where the off-axis pixels contain emission, the dark was obtained from a deep blank-sky observation obtained on the same observing day. We compared the two methods and found no differences, but the noise was smaller using the smooth off axis pixels, in the case of the compact targets.

For all extended sources (Arp 299, ESO 173-G015, MCG+12–02–001, Mrk 331, NGC 1068, NGC 1365, NGC 2146, NGC 3256, NGC 5135, and NGC 7771), an aperture correction is necessary to compensate for the wavelength dependent beam size (Makiwa et al. 2013). We defined a source as extended using LABOCA or SCUBA 350 or 450 μm (respectively) maps with 8'' resolution. We convolved the 8'' resolution maps with the SPIRE FTS resolution, and if the galaxy was more extended than the smallest SPIRE beam size, we defined it as extended. In order to correct for the extended nature of these sources, we employ HIPE's semiExtendedCorrector tool (SECT). This tool "derives" an intrinsic source size by iterating over different source sizes until it finds one that provides a good match in the overlap range of the two observing bands near 1000 GHz, and is further discussed in Wu et al. (2013). We set the Gaussian reference beam to 42'', the largest SPIRE beamsize. The beamsize corrected flux values for the 10 extended sources are listed in Table 2, along with the compact sources. We note that the error in the extended source flux correction could be significant due to the assumptions that the high-J CO transitions are distributed in the same way as the low-J CO lines. If high-J-CO transitions are only coming from a centralized compact region, we are overestimating their flux with our beam correction method. For this reason, we apply an additional 15% error to the extended galaxies. There are three targets in the sample that have multiple pointings; Arp 299, NGC 1365, and NGC 2146. In the case of Arp 299, we use only the pointing for Arp 299 A. For NGC 1365 we take the average of the northeast and southwest pointings. This is done since the northeast and southwest pointings have approximately a 50% overlap in field of view at the center of the galaxy. This overlap region is the center of the PACS observations, so for comparison, it is best to average the northeast and southwest pointings. For NGC 2146, we use only the nuclear pointing.

CO and [C i] line fluxes were extracted using version 1.92 of FTFitter,²⁶ a program specifically created to extract line fluxes from Fourier transform spectrographs, and are listed in Table 2. This is an interactive data language based graphical user interface, that allows the user to fit lines, choose line profiles, fix any line parameter, and extract the flux. We define a third order polynomial baseline to fit the continuum for the SLW and SSW separately and derive the integrated line intensities from baseline subtracted spectra with a simultaneous line fit of all CO, [C i], [N ii] and other bright lines in the spectrum. We use a Gaussian line profile, which is based on the assumed intrinsic line shape with a width derived from spectrally resolved CO 1–0, convolved with the instrumental line shape, which is a sinc profile. We adopt an error of 16% for the non-extended galaxy fluxes, which encompasses our dominant sources of error, 10% for the flux extraction and baseline definition, and 6% for the absolute calibration uncertainty for staring-mode SPIRE FTS observations (Swinyard et al. 2014). For the case of extended sources, we adopt as already mentioned, an additional 15% error from the beam size corrections, resulting in a total error 30% for the 10 extended sources.

We have obtained observations of the [O i] 63 μm ([O i]₆₃), [O i] 145 μm ([O i]₁₄₅), and [C ii] 157 μm emission lines with the Integral Field Spectrometer of the PACS (Poglitsch et al. 2010) on board the Herschel Space Observatory for every object in the HerCULES sample. The data presented here were obtained as part of the Herschel program KPOT_pvanderw_1 (PI: P. van der Werf), complemented by observations from other programs. The observations and program IDs of the [C ii] and [O i] lines are listed in Table 1.

The data were downloaded from the Herschel Science Archive and processed using HIPE v11.0. Standard processing steps including timeline deglitching, application of the Relative Spectral Response Function and detector flat fielding, and subtraction of the on and off chop positions, gridding along the spectral axis, and combination of the nod positions. With the exception of Arp 299, the objects are all centered on the 9 4 central spaxel of the 5 × 5 PACS array, observed in staring mode. The fluxes are extracted from the central spaxel, using the extractSpaxelSpectrum routine, and referenced to a point source. We use the pointSourceLossCorrection routine to capture any additional flux that may not be captured in the central spaxel. Finally, version 3.10 of SPLAT as part of the STARLINK software package (http://star-www.dur.ac.uk/~pdraper/splat/splat-vo/) was used to subtract the baseline from each observation, and isolate the desired lines, in the case of PACS range spectroscopy. The reduction steps were the same for both the PACS range and line spectroscopy, two different observing modes of PACS.

Arp 299 was observed in the mapping mode and reduced using the standard pipeline reduction. The integrated flux for Arp 299 A (presented in this paper) was calculated by summing the flux within a 25'' aperture centered on Arp 299 A SPIRE pointing.

In order to extract the line parameters from the PACS observations, we first integrate over the baseline subtracted spectrum and then we fit a Gaussian profile to the baseline subtracted flux. In some sources, the [O i]₆₃ line shows a double-peaked profile, where the flux at the central wavelength is diminished, which could indicate Keplerian rotation. However, if this were the case, then we would expect a similar profile in the [O i] 145 μm line and possibly the other fine structure lines as well, which is not seen. The spectral resolution of PACS at 145 μm is more than sufficient to resolve the ∼0.2 μm separation between the two peaks in the [O i]₆₃ profile (Figure 1). Therefore, we conclude that this double-peaked profile is due to absorption in the center of the profile by colder foreground gas. We note that [O i]₆₃ absorption is due to O in the ground state while absorption at 145 μm requires O to be at a state having an energy of 226 K above the ground state. Therefore, in cool or moderate density gas the [O i]₁₄₅ line will not show an absorption feature, even if the [O i]₆₃ line does. This same effect has been noted in Arp 220, which shows the [O i]₆₃ in full absorption (González-Alfonso et al. 2012). In the case of NGC 4418 and Zw 049.057, the [O i]₆₃ line has an inverse P Cygni profile, suggesting that the absorbing foreground gas is flowing into the nuclear region. The three example galaxies for which the spectra are shown in Figures 1–3 display increasing absorption of the [O i]₆₃ μm line. In NGC 7552, a face-on starburst galaxy, the profile remains Gaussian, while in Mrk 331, a late-stage merger, there is a strong dip in the middle of the profile. IRAS F17207−0014 is known for being one of the coolest ULIRGs, here absorption dominates the [O i]₆₃ emission. For the [O i]₆₃ profiles that show an absorption feature, we fit the Gaussian only to the wings of the emission profile and state the flux in parentheses. The Gaussian-fit flux is only valid if the true line profile is Gaussian. We suggest this is a more robust estimate of the true integrated flux of the [O i]₆₃ line emerging from the warm nuclear region, since in many cases, the absorption dominates the profile. We note that using a Gaussian profile to extrapolate the line flux requires an assumption that the location of the line center (a free parameter) is in the middle of the profile, which may not be accurate, especially in the case of IRAS F17207−0014, or any other galaxies with an asymmetric line profile. We have tested the relations presented in the rest of this paper with both the integrated flux and the Gaussian fit, and find it does not strongly affect the results. Both the observed line fluxes and the Gaussian-fit line fluxes, stated in the parentheses, are also presented in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All SPIRE CO and [C i] line fluxes are listed in Table 2. We present three examples of galaxy spectra obtained with SPIRE in the top panels of Figures 1–3 for NGC 7552, Mrk 331, and IRAS F17207−0014, respectively. It is important to note that the baseline ripple seen in the SPIRE FTS spectra is due to both the sinc profile of the strong CO transitions and the noise. Since in this paper we only discuss the neutral gas cooling, we do not present fluxes of [N ii] (which originates in ionized gas) or the molecular lines other than CO, which are irrelevant to the total neutral gas cooling. A comprehensive set of fluxes will be presented in P. van der Werf et al. (in preparation). In addition, the HerMES team is planning to publish a formal data paper using HIPE v12.0 with a detailed error analysis in the near future.

In the bottom row of the spectra in Figures 1–3, the PACS line profiles of the three sample galaxies are presented (NGC 7552, Mrk 331, IRAS F17207−0014).

In Figure 4, we present the CO ladders of the full HerCULES sample. We have collected the available ground based observations of CO J = 1–0, 2–1, and 3–2, whose fluxes and references are listed in Table 3. Where necessary we have converted the ground-based measurements to the cosmology adopted here (Section 1). In order to compare these CO ladders directly, we have normalized the ladders by the integrated CO flux summed from J = 4–3 through J = 13–12, to focus on the relative behavior of the higher-J transitions since we do not have CO J = 1–0 data for all sources. For galaxies with non detections for any CO transitions, we use the linearly interpolated value. The CO ladders are separated into three classes based on the parameter α, where:

$$\begin{equation} \alpha =\frac{L_{{\rm CO}_J=11-10}+L_{{\rm CO}_J=12-11}+L_{{\rm CO}_J=13-12}}{L_{{\rm CO}_J=5-4}+L_{{\rm CO}_J=6-5}+L_{{\rm CO}_J=7-6}}. \end{equation} \tag{ 1 }$$

and $L_{{\rm CO}}=4\pi D_L^2 F_{{\rm CO}}$, with F_CO in (W m⁻²) (Table 2) and luminosity distance, $D_L^2$ in (m) (listed in Table 1). Here we use three transitions of both the mid- and high-J CO transitions to help prevent noise or a non-detection of one of these lines from dominating α. We define the three classes as:

• 1.  
  Class I: α < 0.33.
• 2.  
  Class II: 0.33 < α < 0.66.
• 3.  
  Class III: α > 0.66.

Zoom In Zoom Out Reset image size

The definition of the classes is quantitatively arbitrary, but chosen to reflect similarities in the spectral line energy distributions, which is illustrated in Figure 4. In the case that we do not have any observations of the low-J transitions, we do not plot any low-J fluxes. The parameter α is based on the ratio of three high-J CO lines to three mid-J CO lines, which essentially defines the drop-off slope of the CO ladder from J = 5–4. Thus, the steepest drop-offs are in Class I, while the flattest ladders are in Class III. Class II consists of objects that peak around J = 6–5, but do not fall off as steeply as those of Class I. Our three example galaxies were selected to fit into these categories, with NGC 7552 as a Class I, Mrk 331 as a Class II, and IRAS F170207−0014 as a Class III object.

We note that the CO ladders for many Class II and III (the excited classes) objects have been published. In the case of all of these sources, heating mechanisms besides UV heating are required to explain the high-J CO emission, when also considering additional constraints. Arp 220 (Rangwala et al. 2011), Arp 299 (Rosenberg et al. 2014b), NGC 253 (Rosenberg et al. 2014a), and NGC 6240 (Meijerink et al. 2013) require mechanical heating to reproduce the high-J CO lines, while Mrk 231 (van der Werf et al. 2010) and NGC 1068 (Spinoglio et al. 2012) require X-rays to directly heat the gas in order to reproduce the observed molecular emission. This trend suggests that when dealing with highly excited CO ladders, such as in Class II and especially Class III objects, there is an additional heating mechanism necessary to explain the observed molecular emission. We will explore this issue for the full sample in the next section.

In the local universe, the ratio of IRAS 60/100 μm flux densities correlates with infrared luminosity (e.g., Chapman et al. 2003). It is of interest to determine whether the IRAS 60/100 μm ratio (a proxy for dust temperature) or L_FIR (dust luminosity) correlates better with the degree of CO excitation, as parameterized by α. Since α represents a proxy for the slope of the CO ladder above J = 5, it traces the relative brightness of high-J lines in comparison to mid-J lines, allowing a rough estimate of overall CO excitation. When α is small, the CO excitation is low and the CO SLED is highly peaked and when α is large, the CO SLED is flattened and the excitation is high, indicating significant emission by warm and dense molecular gas. We compare the S_60 μm/S_100 μm ratio from Sanders et al. (2003) and the L_FIR to the molecular gas excitation (α). In the left panel of Figure 5 the excitation (α) is plotted as a function of S_60 μm/S_100 μm, where each square point is a galaxy in our sample. The best-fit power law is shown with a red dashed line. In the right panel, a similar correlation between L_FIR and α is shown, with the best-fit power law plotted with a red dashed line. We see that although both panels show a positive trend, the correlation found with the S_60 μm/S_100 μm ratio is tighter than that seen with the L_FIR. The molecular gas excitation to infrared color relation has a correlation coefficient of r = 0.81, while the excitation to L_FIR relation has a correlation coefficient of r = 0.52. Although the excitation to infrared color relationship is significantly correlated, we find three outliers, based on the largest Euclidean distance and shown in red in the left panel of Figure 5), Arp 299, NGC 1614, and NGC 2623. We also select the three farthest outliers, also based on the Euclidean distance, in the α to L_FIR relation (Figure 5, Arp 299, NGC 4418, and IRAS F18293−3413), also plotted as red points in the right panel. We use the traditional definition of Euclidean distance: ||u − v||². To test the strength of the correlations, we refit a power law excluding the three outliers in each plot, which are marked in red. We remove the outliers to test if the correlation between the CO excitation and the IRAS colors is still higher than that of the L_FIR, or if it is just the outliers affecting the coefficient. The new best fit is plotted as the blue solid line in Figure 5. Although the exclusion of these points does not result in a significant change in the best fit in either case, it does improve the correlation coefficients, resulting in a correlation coefficient of r = 0.88 for the molecular gas excitation to infrared color correlation, and r = 0.70 for the molecular gas excitation to L_FIR relation. Physically, this suggests that the presence of warm, dense molecular gas, is correlated with the presence of warm dust. It is, however, important to note that, once removing the outliers, the correlation between α and L_FIR is also significant, with a 1.6% probability of this relation being spurious before removing the outliers. Further, since the IRAS 60/100 μm ratio is shown to correlate with L_FIR, these two quantities are likely related by underlying variables, making this correlation difficult to interpret. We note however that Lu et al. (2014) compared specific CO line transitions normalized by FIR and the IRAS 60/100 μm flux ratio. They find that as the CO gas becomes warmer, the 60/100 μm ratio also increases, which is in good agreement with our results.

Zoom In Zoom Out Reset image size

We can calculate the neutral gas cooling budget in each galaxy by summing the luminosities of the [O i], [C i], [Si ii], [C ii], and CO lines, since these are the main neutral gas coolants in the mid- and far-infrared regime. We take the [Si ii] fluxes from Inami et al. (2013), which were observed with the Spitzer IRS instrument in the long wavelength, high resolution mode. We note that the slit size for Spitzer IRS-LH is 111 × 223, thus galaxies that are more extended than 111 will be missing some [Si ii] flux. There are five affected galaxies, Arp 193, Arp 299, IRAS 13242−5713, NGC 1365, and NGC 2146. In Figure 6, we present the percentage of cooling contributed by each emission line as a function of L_FIR. The percentage cooling for each species is calculated by comparing the luminosity of a species to the total summed luminosity of the [O i], [C i], [Si ii], [C ii], and CO lines. We exclude three galaxies based on their [O i]₆₃ profiles, Arp 220 that is fully in absorption, and NGC 4418 and Zw 049.057, which are both heavily absorbed and show inverse P-Cygni profiles (González-Alfonso et al. 2012). For each galaxy, the percentage of cooling contributed by [C i]₆₀₉ and [C i]₃₇₀ is plotted as a yellow circle, CO J = 4–3 through J = 13–12 in red, [Si ii] in blue, and the combined cooling of [O i] 63 μm, [O i] 145 μm, and [C ii] is plotted in green. In the bottom panel of Figure 6, we separate the cooling contributions of [O i] and [C ii]. For the five galaxies that are larger than 111, we expect a higher [Si ii] contribution than shown in Figure 6. The solid lines show the mean percent-cooling for each emitting species. For example, the neutral atomic carbon is responsible for no more than 2% of the total cooling with an average cooling contribution of 1.5%, while CO contributes a mean of 10.8% and [Si ii] contributes 24.2%. There are two galaxies with exceedingly high CO cooling percentages, namely IRAS F05189−2524 and Mrk 231, the two strongest AGN in the sample. Very high CO cooling percentages have also been noted in the massive Galactic star forming region W3 with 32% total gas cooling (Kramer et al. 2004), and even higher percentages in DR21 (Jakob et al. 2007; White et al. 2010). In both cases, the high CO percentage is attributed to self absorption of the [O i]₆₃ line, yet we calculate our fluxes both with the Gaussian and observed [O i]₆₃ line fluxes and see little change. Therefore, for IRAS F05189−2524 and Mrk 231, we believe the high CO cooling percentage is a true effect, and not one dependent on the [O i]₆₃ absorption. The most efficient coolants are [O i] and [C ii], which together provide a mean of 63.7% of the total gas cooling budget. Separating their cooling contributions, [C ii] provides a mean cooling percentage of 33.6% and [O i] cools a mean of 30.1% of the gas. The mean cooling percentages and their standard deviations are shown in Table 4. Inspection of Figure 6 shows that the outliers are randomly distributed, and there is no clear trend of outliers as a function of L_FIR.

Zoom In Zoom Out Reset image size

Table 4. Mean and Standard Deviation of Percent Cooling Contribution

Line	Mean	Std. Dev.
$[{\rm C\,\scriptsize{II}}] +[{\rm O\,\scriptsize{I}}]_{63+145}$	63.7	14.3
$[{\rm C\,\scriptsize{II}}]$	33.6	9.1
$[\rm O\,\scriptsize{I}]_{63+145}$	30.1	11.8
$[\rm Si\,\scriptsize{II}]$	24.2	9.9
$\sum _{j=4}^{13}\rm CO_{{j}}$	10.8	10.0
$\sum _{j=1}^{2}\rm [{\rm C\,\scriptsize{I}}]_{j}$	1.5	0.9

Download table as:  ASCIITypeset image

It is interesting to note how constant each cooling range is as a function of L_FIR. It is a well-known phenomenon that as far infrared luminosity increases, an apparent [C ii]/FIR ratio decreases, which is what we call the [C ii] deficit. This is observed in various environments in the local universe (Malhotra et al. 2001; Luhman et al. 2003; Díaz-Santos et al. 2013). In addition, Graciá-Carpio et al. (2011) find a fine-structure line deficit in [N ii] and [O i] as well, such that both the [N ii]/FIR and [O i]/FIR ratios decrease as a function of far infrared luminosity. One corollary to this fine-structure line deficit is that other (molecular) coolants could become more efficient at higher L_FIR. However, as shown in Figure 6, the relative efficiencies of each coolant remain mostly constant. In order to better understand this phenomenon, we plot the [C ii] deficit ([C ii] to FIR ratio) as a function of L_FIR in the top left corner of Figure 7. For comparison, we plot the same ratios but for the [O i]_{63 + 145}, [Si ii], CO, and [C i] lines, where the CO lines encompass all transitions 4 ⩽ J_upp ⩽ 13 and [C i] is the sum of both [C i]₆₀₉ and [C i]₃₇₀. We see a trend that is consistent with a deficit in [C ii], agreeing with the observations from Díaz-Santos et al. (2013). Since our sample spans a smaller luminosity range than Díaz-Santos et al. (2013), our deficit only spans a factor of two in the [C ii]/FIR ratio, while in the GOALS sample, it spans an order of magnitude. Thus, we plot the deficit in linear space and see clearly that the only points with a high [C ii]/FIR ratio are those of lower L_FIR. We see a similar trend in [Si ii] and [C i], where the latter has also been observed from ground based observatories (Gerin & Phillips 1998, 2000). The [O i] shows a tentative line deficit, which becomes more obvious when we include galaxies of higher IR luminosities (Graciá-Carpio et al. 2011). Although we find evidence consistent with line deficits in [C ii], [O i], [Si ii], and [C i], there is no evidence for a deficit in CO, which shows a flat distribution over L_FIR, further strengthening the results from Lu et al. (2014). In all panels, NGC 6240 is an extreme outlier, as was also noted in Lu et al. (2014).

Zoom In Zoom Out Reset image size

The difference between the molecular and fine structure emission can be understood in terms of heating mechanisms. The fine structure lines are originating from regions that are heavily affected by UV photons, at the edges of PDRs. As shown in the PDR models from Kaufman et al. (1999), as the radiation field and density increase, the fine structure line emission is expected to weaken compared to the far infrared flux. However, this does not apply for the molecular gas (CO), where we see no line deficit. This result may not be surprising since CO traces the molecular gas deeper in molecular clouds, where the UV field is significantly attenuated. We can test the gas versus dust cooling efficiency by plotting the ratio of the total gas cooling to the far infrared flux as a function of L_FIR, which is shown in the bottom right panel of Figure 7. Here, it is clear that the trend is decreasing in a very similar manner to that of the fine structure lines, which is expected since [C ii], [O i], and [Si ii] dominate the cooling. Therefore, we can say that the fine structure line deficit is actually a gas cooling deficit in comparison to the FIR flux that is a result of UV heating becoming a more efficient coolant of dust in exceedingly extreme environments. It is critical to note that the fine structure line fluxes are not decreasing in absolute flux, but their relative contribution in comparison to the warm dust (measured with the FIR flux) is decreasing, due to the warm dust becoming increasingly efficient at a faster rate than the gas coolants. This can also be understood as dust being heated more efficiently, such that the fine structure lines do not actually show a deficit, which agrees with the results of Díaz-Santos et al. (2013), where they find that grains are "stealing" photons from the gas.

However, if the neutral gas cooling efficiency is decreasing, and we observe a tentative line deficit in all species except for CO, then we would expect the percentage of CO cooling to increase slightly as a function of L_FIR, to compensate for the decreasing fine structure line cooling efficiency. We do see a slight increase in the percentage of CO cooling, but the trend is tentative and within the general scatter of the other CO cooling percentages. In addition, the galaxies with high L_FIR have flatter CO ladders, meaning there is non-negligible flux in high-J (J >13) transitions, that we do not account for in our cooling budget. This missed high-J flux would increase the CO cooling percentage, and could make the CO cooling budget increase for higher luminosity sources. Since we only have four galaxies in our sample within the range that we would expect to see an elevated CO cooling percentage, we would need to increase our sample size in the high luminosity regime to determine if this trend is real.

These results have several implications. First, it is unlikely that the [C ii] line deficit could be caused by high dust opacities at 158 μm, since the line deficit is also observed in lines of much longer wavelengths, notably the [C i] line. Similarly, the line deficits are not likely caused by the line being optically thick, since the deficit is also observed in the optically thin [C i] lines, which agrees with the results of Díaz-Santos et al. (2013). A likely explanation is star formation dominated by ultracompact H ii regions, which would suppress all fine-structure lines due to an increase of dust competition for UV photons. This is further supported by the fact that the total gas cooling ([C ii]+[O i]+CO+[C i]) decreases as a function of L_FIR, meaning that dust becomes an even more efficient coolant as the infrared luminosity increases.

The lack of a strong CO-line deficit shows that the bulk of the molecular gas heating is not affected by the mechanism suppressing the fine structure lines. This is interesting, since in a UV-photon heated environment, suppressing the UV field implies a reduced heating rate, and therefore also a lack of warm molecular gas might be expected. We also note that the integrated CO luminosity, which represents tens of percent of the total gas cooling, is much more than what is predicted by any pure PDR model, which give a CO cooling fraction of at most a few percent (3%–5%) of the total gas cooling (e.g., Meijerink et al. 2011). This result suggests that the CO luminosity may be powered by a different heating mechanism, which does not lead to dissociation or ionization.

The CO molecule can be heated indirectly through (1) the photoelectric effect by ultraviolet (UV) photons, (2) by fast electrons from directly ionized H and H₂ by X-rays, cosmic rays (Meijerink & Spaans 2005), or (3) mechanical processes, which includes shocks and turbulence. X-rays heat gas by ionizing H and H₂ directly, and these fast electrons then thermalize the molecular gas with an efficiency of 10%. UV photons ionize polycyclic aromatic hydrocarbons (PAHs) and dust grains, and the resulting free electrons heat the molecular gas with a net efficiency of 1%–3%. In addition, the chemistry in an XDR is driven by X-ray photons instead of FUV photons that are able to penetrate deeper into the cloud without efficiently heating the dust at the same time. These X-rays are mostly produced by AGNs or in areas of extreme massive star formation. Cosmic rays can also heat the gas by penetrating into cloud centers, similarly to X-rays, and are typically produced by supernovae. Mechanical heating is another efficient source of gas heating. This is commonly attributed to turbulence in the interstellar medium (ISM), which may be driven by supernovae, strong stellar winds, jets, galaxy mergers, cloud–cloud shocks, shear in the gaseous disk, or outflows.

To investigate the main mechanism heating the molecular gas, we use another diagnostic molecule, namely PAHs. PAHs are carbonaceous, nanometer-sized macromolecules that contain 50–100 carbon atoms with an abundance of 10⁻⁷ per hydrogen atom (Tielens 2008). The absorption of one far-UV photon is enough to heat the PAH molecule to a high temperature and will cause this molecule to emit in the characteristic bands at 3.3, 6.2, 7.7, 8.6, and 11.2 μm (Tielens 2008, and references therein). Because PAHs are only fluorescently excited, and are easily destroyed by more energetic radiation, they are ideal tracers of UV heating. Thus, we can use the equivalent width (EW) of the 6.2 μm feature from the Spitzer Space Telescope (Stierwalt et al. 2013), as a proxy for the UV energy density. In Figure 8, we compare the PAH EW to the percentage of the total gas cooling done by CO, as calculated in Section 4.2. Class I objects have the steepest decreasing CO SLED, high PAH EWs, and low percentages of CO cooling. On the other hand, Class II and Class III objects with high percentages of CO cooling have low PAH EWs. We note that the two objects with very low PAH EWs are the most AGN-dominated objects in our sample (Mrk 231 and IRAS F05189−2524), where PAH destruction by X-rays and nuclear hot dust combine to strongly lower the PAH EW. The high CO cooling fractions of these objects can likewise be attributed to energy input by the AGN (through X-ray heating or mechanical energy input from an AGN-driven outflow). There are two galaxies that are Class I objects with lower PAH EWs, and those are IC 1623 and NGC 7469. NGC 7469 is a starbursting galaxy with a Seyfert 1 nucleus and IC 1623 is a late-stage merger with a starburst nucleus. NGC 7469 has a central starburst ring along with its AGN nucleus (e.g., Davies et al. 2007). This suggests that the Class I type CO emission is coming from the starburst ring that is encircling, but not directly affected by, the AGN. These trends reinforce the idea that for objects with a high CO cooling fraction, the CO is efficiently excited by something besides UV photons.

Zoom In Zoom Out Reset image size

We now compare our α parameter with the CO J = 1–0 linewidth for the sources in our sample (Figure 9). We correct the linewidth for inclination such that all galaxies are effectively turned edge on using the K-band axis ratio from Two Micron All Sky Survey, with the exception of Arp 220 and NGC 6240. For Arp 220 we used the inclination derived by Scoville et al. (1997) using arcsecond imaging of CO J = 1–0 and for NGC 6240 we used the inclination derived by Engel et al. (2010) using subarcsecond near-infrared imaging. The range of Class I objects is highlighted with a red, Class II with a green, and Class III with a blue background.

Zoom In Zoom Out Reset image size

Merging and interaction, molecular outflows, and random motions may all play a role in increasing the linewidths, but are not expected to dominate the linewidths of our target sources. The linewidth of CO J = 1–0 is dominated by rotation, which is determined by the mass of the central regions of the galaxy. Therefore, galaxies with high excitation (α) and high linewidths are more massive. However, the speed of rotation is also a proxy for the mechanical energy reservoir available in the galaxy nucleus since rotation leads to processes such as sheering and turbulence. If the fraction of mechanical energy that is converted into heating the molecular gas is constant, then high linewidth galaxies in Figure 9 would have more mechanical energy available, and would be responsible for heating more of the molecular gas. However, establishing this result would require a detailed study of the velocity fields of the molecular gas in all of our targets, now possible with ALMA. Since such data is not available, thus we suggest that the linewidth provides an estimate of the available reservoir of mechanical energy, but we remain agnostic as to the extent that this reservoir is actually tapped.

Inspecting Figure 9 with these considerations in mind, the fact that low-excitation galaxies all have the smallest linewidths, while high-excitation galaxies have a range of linewidths has an attractive physical interpretation. In high excitation galaxies with high linewidths, mechanical energy input is a viable source of excitation, while in high excitation galaxies with low linewidths, radiative energy input may be more important. Indeed, in Figure 10, the lowest linewidths among high excitation galaxies are found in Mrk 231 and NGC 4418 (an exposed and an obscured AGN, respectively).

Zoom In Zoom Out Reset image size

Class I and II objects with low linewidths are likely dominated by UV heating, since their CO ladders turn around somewhere before J = 7–6. For comparison, we have highlighted the position of our three example galaxies, NGC 7552 (Class I), Mrk 331 (Class II), and IRAS F17207−0014 (Class III), as well as some of the more famous galaxies in our sample. The major merger galaxies NGC 6240 and Arp 220 have relatively high α values and linewidths around 650 km s⁻¹. On the other hand, the strongest AGN in our sample, Mrk 231, lies in the low linewidth region of Class III, suggesting that its gas is radiatively heated, likely by X-rays from the AGN. This confirms the results of van der Werf et al. (2010), where they find the high-J CO excitation consistent with an XDR. The other two galaxies with AGN contribution are NGC 1068 and IRAS F05189−2524, both of which are Class II objects that lie in the low linewidth region, suggesting that they are also radiatively excited. However, since both of these objects lie in Class II (albeit on the border between Class II and III), it is unclear whether their excitation is from UV, higher energy photons, or cosmic rays. Both NGC 1068 and IRAS F05189−2524 have also been studied in detail by Spinoglio et al. (2012) and Pereira-Santaella et al. (2014), respectively. For NGC 1068, Spinoglio et al. (2012) and Garcia-Burillo et al. (2014) find that indeed, the excitation is due to either XDRs (in the circumnuclear disk) or PDRs (in the star-forming ring), which agree with our results. In IRAS F05189−2524, Pereira-Santaella et al. (2014) find that there is a large contribution from mechanical heating in this source, as traced by a shallow H₂ temperature distribution, yet they cannot rule out the AGN as an important heating source for the molecular gas.

In order to examine this method vis-a-vis the role of AGNs and to determine additional heating sources, and to check if there are any underlying biases in the J = 1–0 linewidth, we can compare the inclination corrected linewidths with the percentage AGN contribution for each galaxy. In the case of very disturbed major mergers, such as Arp 220, Arp 299, IRAS 17208−0014, and NGC 6240, the inclination is difficult to determine, and thus the corrected linewidth can be over- or underestimated. The uncorrected linewidths are shown in Table 1 for comparison. The AGN contribution can be estimated using both the 15/30 μm flux density ratio (f15/f30) and the [Ne v]/[Ne ii] ratio along with the prescriptions from Veilleux et al. (2009). We take the f15/f30 ratio from Stierwalt et al. (2013) and the [Ne v]/[Ne ii] ratio from Inami et al. (2013). We average the results of the AGN contribution from both methods and find average AGN contributions ranging from to 0%–95% of the bolometric luminosity. For the high excitation sources (Class III) we use the average AGN contribution in combination with the inclination corrected linewidth to separate mechanical heating from AGN heating. In Figure 10, we plot the AGN contribution on the x-axis and the inclination corrected CO J = 1–0 linewidth along the y-axis. This scatter plot shows that for our high-excitation galaxies, very high linewidths are not associated with high AGN contributions and conversely the galaxies with high AGN contributions do not display high linewidths. We also plot the Class II galaxies as green points (the template galaxies are marked) for comparison. Most Class II galaxies have low AGN contributions, while they have a large range of linewidths.

For example, we can compare Arp 220 and Mrk 231. Both galaxies have high α values, but Arp 220 has a high linewidth and a low AGN contribution, while Mrk 231 has a low linewidth and a high AGN contribution. Comparing NGC 6240 and Arp 220 reveals that although both have high linewidths and α values, NGC 6240 also has a high AGN contribution. This suggests that both the AGN and mechanical processes are contributing to heating the gas in NGC 6240.

It is difficult to conclude anything definitive about objects that show average or typical values of α or linewidths. It is also important to note that although there is a non-negligible contribution of heating from mechanical processes or the AGN in Class III galaxies, the gas is still mostly heated through UV heating processes. We caution the use of any two of these diagnostics alone, since for example, a Class III object with a low linewidth may still be mechanically heated by small-scale turbulence that would not produce an observable global line broadening effect. In this case, the AGN contribution would be low, but the α would be high, discounting AGN heating. Using all three parameters simultaneously allows for a qualitative estimate of which additional processes are exciting the warmest molecular gas. In order to fully understand the excitation mechanisms and physical parameters of the molecular gas, an additional detailed modeling of the ¹²CO, ¹³CO, and dense gas tracers (HCN, HNC, HCO⁺, etc.) is required (e.g., Rosenberg et al. 2014a).